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ABSTRACT 


This  report  develops  a  nonlinear  model  which  can  be 
used  to  predict  combustion  instability  zones  in  liquid 
rocket  engines.  The  nonlinear  model  is  developed  by  combin¬ 
ing  a  non  linear  instability  model  with  a  steady -state  vapori¬ 
zation  model.  Such  an  analysis  determines  the  zones  of  an 
engine  in  which  a  tangential  mode  of  high  frequency  insta¬ 
bility  is  most  easily  initiated.  A  rocket  engine  can  be 
analyzed  by  incrementally  dividing  the  combustion  chamber 
in,to  annular  nodes  in  the  r  and  z  directions.  Steady-state 
properties  at  each  annular  node  or  position  in  the  chamber 
are  computed  from  the  steady-state  vaporization  computer 
program.  The  steady-state  program  is  capable  of  computing 
combustion  profiles  in  thermally  unstable  propellants  of 
the  monomethy lhydrazine/nitrogen  tetroxide  type.  This  model 
describes  droplet  vaporization  with  vapor  phase  decomposition. 
Using  the  computed  steady-state  properties  and  the  stability 
limit  curves  from  the  instability  computer  program,  stability 
at  each  node  is  determined.  This  process  is  repeated  for 
each  node  to  determine  a  stability  map  of  the  entire  engine. 
Thus  stability  can  be  related  to  hardware  design  parameters, 
thereby  enabling  the  influences  of  system  design  and  stability 
rating  devices  to  be  determined. 

This  report  has  been  divided  into  four  main  parts.  Part 
one  deals  with  the  application  of  the  steady-state  and  insta¬ 
bility  program  to  determine  stability  zones  in  a  rocket 
engine.  Part  two  covers  the  details  of  the  steady-state 
model  for  monopropellant  type  fuels  such  as  monomethyl- 
hydrazine.  Part  three  deals  with  the  nonlinear  instability 
model  used  to  generate  the  stability  limit  curves  and 
finally,  part  four  contains  the  details  of  the  steady-state 
and  instability  programs  and  computer  listings. 
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I.  INTRODUCTION 


The  problem  of  unstable  combustion  in  rockets  represents 
one  of  the  most  serious  obstacles  to  reliable  performance. 

Despite  the  research  which  has  been  performed  in  an  attempt 
to  understand  the  mechanism  of  combustion  instability,  there 
has  not  been  a  satisfactory  model  relating  stability  to  pro¬ 
pellant  properties  and  combustion  chamber  design.  The  research 
developed  herein  is  directed  at  obtaining  an  understanding  of 
the  behavior  of  propellant  sprays  in  a  liquid  propellant  engine 
under  oscillating  flow  conditions  as  a  means  of  relating  com¬ 
bustion  stability  to  propellant  properties  and  injector  design. 

A  number  of  different  types  of  instability  are  recognized 
which  differ  in  observed  characteristics  and  which  depend  on 
different  mechanisms.  The  work  discussed  in  this  report  is 
concerned  only  with  liquid  propellant  rocket  engines.  Further, 
of  the  various  types  of  instabilities  observed  in  liquid  pro¬ 
pellant  engines,  this  work  is  pertinent  to  the  type  usually 
defined  as  the  high  frequency  instability.  The  principal 
characteristics  of  such  an  instability,  in  terms  of  the  theore¬ 
tical  analysis,  is  that  the  processes  involved  in  these  insta¬ 
bilities  are  assumed  to  be  isolated  within  the  combustion  chamber 
and  are  not  coupled,  nor  do  they  irteract  with  processes  outside 
the  combustion  chamber. 

The  steady-state  combustion  process  can  be  considered  to 
consist  of  the  following  processes  occurring  in  sequerce  or 
in  parallel:  (1)  atomization;  (2)  vaporization;  (3)  mixing; 
and  (4)  chemical  reaction.  Certain  assumptions  and  approxi¬ 
mations  must  be  made  in  order  to  treat  the  problem  theoreti¬ 
cally  without  introducing  such  a  great  complexity  that  useful 
solutions  cannot  be  obtained.  In  order  to  define  the  theore¬ 
tical  approach  and  objectives  of  the  study  reported  here,  a 
qualitative  discussion  of  the  combustion  process  in  the  com¬ 
bustion  chamber  is  useful. 

For  the  steady-state  combustion  process  as  well  as  the 
unstable  combustion,  it  is  assumed  that  atomization  of  the 
propellant  occurs  very  rapidly  close  to  the  injector  face. 

This  does  not  mean  that  atomization  is  ignored.  On  the  con¬ 
trary,  it  is  assumed  that  the  atomization  process  determines 
the  size,  position  and  velocity  distribution  of  the  drops 
which  participate  in  the  combustion  process.  Changes  in 
conditions  which  affect  atomization  thus  affect  the  droplet 
distribution.  The  description  of  the  atomization  process 
becomes,  since  it  is  omitted  from  the  theory,  an  essential  part 
of  an  experimental  program.  The  experimental  program  determines 
the  relationship  between  the  injector  variables  and  droplet 
distributions.  The  theory  begins  with  a  droplet  distribution. 
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Droplet  vaporization  is  treated  in  the  theory  as  the 
rate  controlling  process.  This  implies  that  the  rates  of 
mixing  and  chemical  reaction  are  extremely  rapid  compared  to 
vaporization  rates.  On  such  a  basis,  vaporization  represents 
a  major  part  of  the  total  process  and  only  a  small  error  is 
introduced  by  neglecting  the  time  for  mixing  and  chemical 
reaction.  The  work  of  Pricm  and  co-workers  (Kef.  1)  indicates 
considerable  justification  for  such  an  assumption  for  many 
propellant  combinations.  If  the  mixing  and  chemical  reaction 
times  are  of  the  same  order  of  magnitude  as  the  vaporization 
times,  or  a  more  extreme  case,  if  the  vaporization  time  is 
short  compared  to  the  time  required  for  the  other  processes, 
the  model  does  not  apply.  The  assumption  of  a  combustion 
process  controlled  by  vaporization  rate  thus  becomes  an 
essential  part  of  the  theoretical  discussion  which  follows. 

The  theory  does  not  apply  to  propellant  which  evaporates  so 
rapidly  that  the  combustion  process  occurs  essentially  in  the 
gas  phase  or,  of  course,  for  propellants  injected  as  gases. 

It  is  important  to  remember,  however,  that  only  one  component 
of  a  bi-propellant  system  need  have  a  slow  vaporization  rate 
for  the  assumption  to  apply.  For  steady-state  combustion  it 
is  assumed  that  a  given  drop  distribution  enters  at  the  injector 
face,  and  that  the  vaporization  process,  including  rapid  mixing 
and  reaction,  occurs  within  the  combustion  chamber. 

The  previous  discussion  describes  a  model  for  steady-state 
combustion.  Presumably,  once  initiated,  such  a  system  should 
continue  burning  in  a  stable  manner.  The  problem  of  combustion 
instability  becomes  involved  when  a  disturbance  produces  a 
change  in  distribution,  pressure,  temperature,  velocity  or 
some  combination  of  these.  While  it  is  possible  to  visualize 
a  combustion  process  which  is  unstable  by  itself  (linear 
instabi lity) , the  theory  described  here  requires  some  initiating 
disturbance  (nonlinear  instability).  Such  a  disturbance  may 
be  a  starting  transient,  a  transient  variation  in  injection, 
or  any  change  in  the  established  steady-state  conditions.  The 
interaction  between  such  a  disturbance  and  the  normal  combus¬ 
tion  process  may  lead  to  a  decrease  or  increase  in  the  intensity 
of  the  disturbance.  A  decrease  in  intensity  would  result  in 
a  return  to  stable  operation,  while  an  increase  in  intensity, 
if  continued,  could  lead  to  a  disturbance  level  which  results 
in  engine  failure. 

An  initial  disturbance  may  build  up  to  intensity  levels 
or  originate  at  an  intensity  level  which  completely  changes 
the  combustion  process.  For  example,  large  amplitude  pressure 
waves  could  shatter  droplets  producing  small  droplets  that 
change  the  model.  The  existence  or  development  of  a  strong 
pressure  wave  could  also  lead  to  detonation  combustion.  The 
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theory  developed  in  this  report  does  not  describe  such  pheno¬ 
mena.  The  work  conducted  by  Dynamic  Science  Corporation  is 
divided  into  two  models:  (1)  linear  model;  and  (2)  nonlinear 
model.  The  linear  model  describes  the  amplification  of  an 
initially  sm  11  disturbance  which  could  lead  to  other  processes. 
In  this  sense  the  linear  model  is  concerned  with  the  initial 
phases  in  the  transition  from  a  small  disturbance  to  instability. 
The  linear  model  is  described  in  detail  in  Reference  2,  and 
will  not  be  repeated  in  this  report.  The  nonlinear  model  is 
concerned  with  the  transition  from  a  finite  disturbance  to  a 
large  amplitude  disturbance. 

Combustion  instability  zones  in  a  liquid  rocket  engine 
arc  determined  in  this  report  by  combining  a  nonlinear  insta¬ 
bility  program  with  a  steady-state  vaporization  program, 
providing  an  analytical  framework  for  determining  the  relation¬ 
ship  of  design  parameters  to  stability.  The  analysis  determines 
the  zones  of  an  engine  in  which  a  tangential  mode  of  high  fre¬ 
quency  instability  is  most  easily  initiated.  The  Dynamic 
Science  Corporation  instability  program  considers  the  nonlinear 
conservation  equations  with  mass  addition  using  a  steady-state 
vaporization  expression  for  the  burning  rate.  From  such  a 
model,  important  nonlinear  phenomena  arc  predicted,  i.e., 

(1)  stability  dependence  on  disturbance  amplitude,  (2)  the 
limiting  amplitude  of  pressure  oscillations,  and  (3)  non- 
sinusoidal  waveforms.  This  model  applies  to  a  one-dimensional 
annulus  of  small  length  (Az)  and  thickness  (Ar)  . 

Applying  the  results  of  this  model,  a  rocket  engine  can 
be  analyzed  by  incrementally  dividing  the  combustion  chamber 
into  annular  nodes  in  the  r  and  z  directions.  Steady-state 
properties  at  each  annular  node  or  position  in  the  chamber  are 
computed  from  the  steady-state  vaporization  computer  program. 
These  steady-state  properties  and  the  stability  limit  curves 
from  the  instability  computer  programs  are  used  to  determine 
the  stability  at  each  node.  This  process  is  repeated  for 
each  node  to  determine  a  stability  map  of  the  entire  engine. 

Thus  stability  can  be  related  to  hardware  design  parameters, 
i.e.,  parameters  related  to  injector  design,  chamber  con¬ 
figuration,  and  propellants,  thereby  enabling  the  influences 
of  system  design  and  stability  rating  devices  to  be  determined. 

The  following  report  has  been  divided  into  four  main  parts. 
Part  one  deals  with  the  application  of  the  steady-state  and  in¬ 
stability  program  to  determine  stability  zones  in  a  rocket 
engine.  Part  two  covers  the  details  of  the  steady-state  model 
for  monopropellant  type  fuels  such  as  monomethy lhydrazine . 

Part  three  deals  with  the  nonlinear  instability  model  used  to 
generate  the  stability  limit  curves  and  finally,  Part  four 
contains  the  details  of  the  steady-state  and  instability  pro¬ 
grams  and  computer  listings. 
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II.  APPLICATION  OP  STEADY-STATE  AND 
COMBUSTION  INSTABILITY  MODEL 


1.  Introduction. 


Dynamic  Science  Cor  oration  has  developed  a  nonlinear 
model  for  determining  the  zones  of  a  liquid  rocket  engine  in 
which  a  tangential  mode  of  high  frequency  instability  is  most 
easily  initiated.  This  model  has  been  related  to  hardware 
design  parameters,  i.e.,  parameters  related  to  injector  design, 
chamber  configuration,  and  propellants,  thereby  enabling  the 
influences  of  system  design  and  stability  rating  devices  to  be 
determined . 

A  method  for  determining  the  zones  of  a  liquid  rocket 
engine  in  which  a  tangential  mode  of  high  frequency  is  most 
easily  initiated  was  developed  by  Beltran  and  Frankel  (Ref.  3). 
This  method  uses  the  Priem-Guentert  (Ref.  4)  nonlinear  model 
in  conjunction  with  the  Priem-Heidmann  (Ref.  1)  propellant 
vaporization  model.  A  rocket  engine  is  analyzed  by  incrementally 
dividing  the  combustion  chamber  into  annular  nodes  in  the  r  and  z 
directions  as  illustrated  in  Figure  1. 

From  a  nonlinear  model  important  nonlinear  phenomena  are 
predicted,  i.e.,  (1)  stability  dependence  on  disturbance  wave 
shape,  amplitude,  type  (velocity  or  pressure),  and  position; 

(2)  the  limiting  amplitude  of  the  unstable  pressure  oscillations; 
and  (3)  the  shape  of  the  unstable  wave  forms.  Such  information 
is  not  only  of  use  from  the  preliminary  design  standpoint,  but 
also  as  an  invaluable  tool  in  understanding  how  rocket  engines 
should  be  disturbed  during  their  development  test  program  to 
determine  the  degree  of  stability  of  an  engine.  The  Dynamic 
Science  Corporation  model  enables  the  engine  designer  to  deter¬ 
mine  the  position  to  introduce  the  disturbance,  a  reasonable 
disturbance  amplitude  criteria,  and  the  most  effective  wave 
profile  (pressure  and  velocity  disturbance).  Since  the  injec¬ 
tor  design  variables  can  oe  related  to  thresnoid  disturoance 
amplitude  parameters  can  be  modified  to  increase  stability  of 
a  given  engine  configuration  or  be  used  in  the  preliminary 
design  of  an  engine. 

While  the  Priem-Guentert  instability  model  (Ref.  4)  solved 
the  nonlinear  conservation  equations  with  a  steady-state 
vaporization  expression  for  the  burning  rate,  the  results  are 
only  valid  for  a  large  droplet-gas  relative  velocity  or  drop¬ 
let  Reynolds  number  based  on  the  speed  of  sound.  The  solution 
of  the  nonlinear  model  (Ref.  4)  used  an  explicit  first  order 
finite  difference  scheme.  Different  integration  techniques 
were  attempted  for  the  solution  of  the  nonlinear  model  in 
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FIGURE  I.  Transtage  Rocket  Engine  Model 
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Reference  5,  Difficulties  were  encountered  in  attempting 
to  reproduce  the  solutions  of  Reference  4.  It  was  detcr- 
minded  that  numerical  instability  was  generating  false 
indications  of  combustion  instability  in  the  nonlinear 
model  for  small  values  of  the  burning  rate  parameter  L. 

Under  NASA  con tract  NAS  3-b7  7  the  Priem-Guentert  results 
(Ref.  4)  were  re  calculated  by  Belt  ran  and  Wright  (Ref.  6) 
using  predictor- corrector  formulae.  Stability  limit 
curves  were  generated  to  determine  the  influences  of  droplet 
Reynolds  number.  Results  from  this  program  will  be  used  for 
the  following  report. 

The  steady-state  combustion  model  being  used  for  this 
study  is  a  s teady- s t at e  ,  adiabatic,  one- dimensions  1  flow  model 
treating  propellant  evaporation  as  the  rate  limiting  step  in 
the  combustion  process.  The  program  treats  storable  propellants 
by  using  a  two-flame  model  to  account  for  propellant  decomposi¬ 
tion.  Such  a  program  is  unique  in  that  it  can  calculate  steady- 
state  combustion  rates  of  hydrazine  type  propellants.  Steady- 
state  properties  at  each  annular  node  or  position  in  the 
chamber  are  computed  from  the  propellant  vaporization  program. 
These  steady-state  properties  and  the  curves  from  the  insta¬ 
bility  program  were  us  d  to  determine  the  stability  of  the 
node.  This  process  is  repeated  for  each  node  to  determine  a 
stability  map  of  the  entire  engine. 

Major  results  of  this  analysis  show  that  the  amplitude 
and  position  of  a  pressure  disturbance  required  to  initiate 
instability  can  be  determined,  thereby  defining  a  sensitive 
zone  (and  the  best  place  to  disturb  the  engine).  This  sensitive 
zone  extends  several  inches  from  the  injector  face  and  occurs 
where  the  average  droplets  are  moving  the  slowest  relative  to 
the  gases. 

Dynamic  Science  Corporation  has  shown  that  an  annular 
combustor  section  will  be  more  stable  as  the  droplet  Reynolds 
number  approaches  zero.  Thus,  a  significant  result  of  this 
work  is  that  there  are,  for  a  vaporization  controlled  combus¬ 
tion  process,  three  parameters  affecting  stability: 

(a)  Burning  rate  parameter  -  L 

(b)  Absolute  value  of  relative  velocity  -  Av 

(c)  Reynolds  number  of  drop  based  on  speed  of  sound  -  Re^ 
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2. 


Nodal  Method. 


Combustion  instability  zones  in  a  liquid  rocket  engine 
are  determined  by  combining  a  nonlinear  instability  model  with 
a  steady-state  vaporization  program.  The  analysis  determines 
the  zones  of  an  engine  in  which  a  tangential  mode  of  high 
frequency  instability  is  most  easily  initiated.  In  addition, 
it  represents  an  a  priori  method  for  determining  the  relation¬ 
ship  of  design  parameters  to  stability. 


The  instability  model  considers  the  nonlinear  conservation 
equations  with  mass  addition  using  a  steady-state  expression 
for  the  burning  rate.  This  model  applies  to  a  one-dimensional 
annulus  of  small  length  (Az)  and  thickness  (Ar)  shown  in 
Figure  1.  Applying  the  results  of  this  model,  a  rocket  engine 
can  be  analyzed  by  incrementally  dividing  the  combustion  chamber 
into  annular  nodes  in  the  r  and  z  directions.  Steady-state 
properties  at  each  annular  node  or  position  in  the  chamber  are 
computed  from  a  steady-state  vaporization  program.  These  steady- 
state  properties  and  the  curves  from  the  instability  model  are 
used  to  determine  the  stability  of  the  node.  This  process  is 
repeated  for  each  node  to  determine  a  stability  map  of  the 
entire  engine. 


For  a  vaporization  controlled  combustion  process,  the 
significant  parameters  affecting  .stability  are  L,  Av,  and  Red 
The  burning  rate  parameter  (L)  is  derived  from  Damkohler's 


similarity  group  based  on  the  speed  of  sound 


The  Reynolds  number  of  the  droplet  based  on  the  speed  of  sound 


is  derived  from  the  droplet  vaporization  equation. 


For  a  steady-state  combustion 
are  functions  of  r,  6,  and  z. 
sidered,  these  properties  are 
may  be  defined  as: 

process, Av,  ,  o0,  r 

If  one-dimensional  f 
only  functions  of  z. 
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Therefore,  since  ■  and  A  are  functions  of  z  only,  L  is  a 
function  of  r  and  z.  Since  annuli  are  being  considered  in 
this  model  it  is  assumed  that  the  contraction  ratio  of  any 
annulus  at  z  is  equal  to  the  contraction  ratio  of  the  chamber 
at  z.  The  Reynolds  number  of  the  droplet  is  a  function  of  z. 

The  steady-state  vaporization  program  was  used  to  compute 
Av(z),  Re^z)  and  m(z).  This  program  assumes  that  vaporization 
is  the  controlling  combustion  process,  which  is  equivalent  to 
assuming  that  mixing  and  reaction  rates  are  fast  compared  to 
vaporization  rates  and  that  reacted  products  are  formed  at  the 
same  rate  as  the  propellants  are  vaporized.  This  condition 
is  satisfied  in  most  high  performance  rocket  engines  where 
the  propellant  is  injected  uniformly  over  the  entire  injector 
face.  Since  it  is  also  assumed  that  the  combustion  processes 
begin  at  the  point  where  droplets  are  formed,  a  length  required 
for  atomization  was  added. 

Example 

To  illustrate  the  nodal  method,  combustion  of  a  mono¬ 
methyl  hydrazine  spray  in  nitrogen  tetroxide  was  considered. 

This  method  can  be  used  on  other  propellant  combinations  as 
well.  Calculations  were  based  on  a  monodispersed  spray  with 
a  drop  radius  of  0.093  inches  (75  microns). 

Typical  transtage  engine  parameters  were  used:  chamber 
diameter,  11. 6S  inches;  chamber  pressure,  100  psia;  chamber 
contraction  ratio,  2.43;  initial  drop  velocity,  1000  inches 
per  second;  and  initial  drop  temperature,  530*R. 

Using  the  outlined  chamber  parameters  Avz,  Re<j,  and  m 
are  computed  as  functions  of  z  and  plotted  in  Figures  2,  3  v 

and  4  respectively.  Using  the  vaporization  program,  only 
the  relative  velocity  in  the  z  direction  is  computed;  however, 
in  the  injection  zone  there  are  large  gradients  and  recircu¬ 
lation  zones  due  to  the  viscous  action  between  the  sprays 
and  combustion  gases,  so  that  velocity  differences  exist 
in  the  r  and  6  directions  as  well. 

A  level  of  turbulence  equivalent  to  of  0.01  is  reasonable 

ao 

in  view  of  the  measurements  made  by  Hersh . (Ref . 7)  If  it  is  assumed 
that  the  drops  are  unaffected  by  turbulent  oscillations,  which 
is  justified  with  the  drop  size  group  considered,  this  velocity 
can  be  added  to  AVZ  thereby  obtaining  the  total  Av.  A  mean 
drop  size  was  selected  to  represent  the  average  Av  since  the 
smallest  droplets  will  follow  the  gas  velocity  and  the  largest 
droplets  will  lag.  The  choice  of  a  mean  drop  size  serves  to 
illustrate  the  method  used;  however,  this  does  not  imply 
necessarily  that  a  mean  drop  is  a  good  representation  of  the 
spray  drop  distribution.  The  method  is  quite  flexible  in 
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FIGURE  4.  Profile  of  Steady-State  Vaporization  Rate 


that  it  can  be  readily  extended  to  treat  distributed 
relative  velocities. 

In  Figure  5,  L  is  plotted  versus  r  and  z.  Since  L  is 
directly  proportional  to  r  (equation  II-2),  it  approaches  zero 
at  the  centerline  of 

pressure  amplitude 
values  of  the  parame 
equilibrium  pressure  amplitude  that  the  wave  will  approach  at 
steady-state  and  the  lower  boundary  represents  the  minimum 
pressure  amplitude  necessary  to  initiate  instability.  These 
curves  were  plotted  on  the  basis  of  Priem's  stability  limit 
curves.  Figures  2-6  were  cross-plotted  to  determine  pulse 
pressure  contour  lines  in  the  r- z  plane  in  Figure  7.  These 
isobaric  lines  describe  the  tangential  threshold  pulse 
required  to  initiate  instability.  It  can  be  seen  that  the 
sensitive  zone  is  from  0.2  to  3.0  inches  from  the  injector 
face.  For  the  case  considered,  this  zone  is  largest  at  the 
outer  radius.  In  general,  the  largest  zone  will  occur  where 
Re^  and  Av  determine  a  point  on  the  locus  of  minima  of  Figure  6. 
Thus",  for  a  larger  diameter  engine  the  maximum  sensitive  zone 
may  occur  at  some  intermediate  radius.  A  plot  of  AP  required 
to  induce  instability  versus  z  at  the  outer  radius  is  shown  in 
Figure  8.  This  plot  shows  that  there  is  an  extremely  sensitive 
region  occurring  when  the  relative  velocity  between  droplets 
and  gas  reverses  direction,  i.e.,  when  Av  is  minimum. 

The  method  presented  can  be  very  useful  in  the  design  of 
stable  rocket  engines.  In  connection  with  engine  design,  it 
is  useful  to  know  the  effect  of  various  design  parameters  on 
the  shape  and  extent  of  the  sensitive  zone.  Although  the 
analysis  presented  here  can  be  extended  to  treat  a  variety 
of  design  conditions,  certain  qualitative  observations  are 
immediately  evident.  The  most  sensitive  region  can  be  moved 
closer  to  the  injector  face  by  (a)  decreasing:  1)  mass-median 
drop  Radius,  2)  injection  velocity,  3)  chamber  contraction 
ratio,  or  (b)  increasing:  1)  droplet  geometric  standard 
deviations,  2)  propellant  injection  temperature,  and  3) 
propellant  volatility.  The  sensitive  region  can  be  moved  away 
from  the  injector  by  opposite  parameter  changes.  Since  it  is 
difficult  to  change  one  of  the  above  paraaete/s  independently 
of  the  others,  it  is  recommended  that  ai:y  proposed  design 
change  be  analyzed  in  detail  by  the  methods  presented  here 
before  any  conclusions  are  drawn. 


^the  engine.  Figure  6  is  a  plot  of  the 

required  to  sustain  a  wave  for  various 
W  AV.  The  upper  boundary  represents  the 
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IS 


PRESSURE 


U 


Distance  froa  Injector  -  Inches 


3. 


Conclusion . 


Although  the  results  presented  here  are  reasonable  in 
light  of  liquid  rocket  engine  testing  experience,  it  must  be 
kept  in  mind  that  the  basis  of  the  model  is  a  one- dimens ional 
description  of  the  transient  processes.  Two-  and  three-di¬ 
mensional  solutions  would  provide  a  description  of  the  multi¬ 
dimensional  nonlinear  wave  interactions  and  permit  a  more 
realistic  consideration  of  the  chamber  boundary  conditions 
which  may  significantly  affect  the  results.^  However,  it  is 
felt  that  the  agreement  noted  between  the  one-dimensional 
model  and  experimental  findings  (Ref.  8  and  9)  indicate 
that  the  one- dimens i oiu  1  approximation  leads  to  meaningful 
results.  Thus,  the  model  is  valuable  from  an  engineering 
standpoint  since  it  relates  design  parameters  directly  to 
stability  determination. 

Results  of  this  study  show  that:  1)  the  amplitude  and 
position  of  a  pressure  disturbance  required  to  initiate 
instability  can  be  determined,  thereby  defining  a  sensitive 
zone;  2)  this  sensitive  zone  extends  several  inches  from  the 
injector  face;  3)  the  most  sensitive  region  in  the  zone  occurs 
where  the  average  droplets  are  moving  the  slowest  relative  to 
the  gases;  and  4)  approximate  baffle  length  can  be  determined, 
i.e.,  baffles  must  be  extended  past  the  instability  zone. 
Dynamic  Science  Corporation  has  shown  that  as  the  droplet 
Reynolds  number  approaches  zero,  a  rocket  engine  becomes 
more  stable.  A  significant  result  of  this  work  is  that  there 
are,  for  a  vaporization  controlled  combustion  process,  three 
parameters  affecting  stability:  1)  burniifc  rate  purs ieter; 

2)  absolute  values  of  the  relative  velocity'  and  3)  Reynolds 
number  of  drop  based  on  speed  of  sound-Re  . 
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III.  STEADY-STATE  DROPLET  COMBUSTION  IN  THE  MONO 
ME THY LHY DRAZ I NE- NITROGEN  TETROXIDE  SYSTEM 


1 .  Steady-State  Vaporization  with  Vapor  Phase  Decomposition. 

■'.lalysis  of  high  frequency  instability  in  liquid  rocket 
engines  requires  a  spatial  knowledge  of  such  parameters  as 
propellant  burning  rate,  relative  velocity  between  gases 
and  propellant  drops,  and  Reynolds  number  of  propellant  drops 
based  on  the  local  speed  of  sound.  This  section  describes  a 
steady-state  combustion  model  which  estimates  these  para¬ 
meters  in  the  combustion  chamber. 

A  model  for  steady-state  droplet  combustion  allowing  for 
vapor  phase  decomposition  of  thermally  unstable  fuels  has 
been  developed  and  will  be  discussed.  In  this  report  the 
model  was  applied  to  the  monomethy lhydrazine/nitrogen  tetroxide 
system.  The  model  presented  enables  the  computation  of  steady- 
state  profiles  within  a  rocket  combustion  chamber.  From  these 
profiles  combustion  instability  parameters  can  be  computed. 
Details  of  the  computer  program  developed  from  the  steady- 
state  model  are  presented  in  Appendix  I.  With  injection 
parameters  specified,  i.e.,  droplet  size,  injection  velocity, 
chamber  pressure,  mixture,  ratio,  and  droplet  temperature,  the 
computer  program  can  compute  droplet  and  gas  histories  along 
the  chamber.  This  information  is  required  for  the  design  of 
combustion  chambers  and  calculation  of  C*  efficiency.  From 
the  results  of  the  computer  program  the  effects  which  any 
of  the  numerous  design  or  physical  parameters  have  upon 
chamber  operation  can  be  determined.  Thus  the  influences 
of  injector  design  parameters  and  physical  properties  on 
combustion  instability  can  be  evaluated. 

This  program  assumes  that  droplet  vaporization  is  the 
rate  controlling  step  in  liquid  rocket  chamber  combustion. 

This  approach  has  proven  successful  in  previous  studies, 

(Ref.  1,  10).  The  vaporization  rate  of  a  droplet  is  deter¬ 
mined  by  transport  phenomena  across  a  boundary  layer  diffu¬ 
sion  mantle.  For  such  fuels  as  MMH  and  NTO,  decomposition 
occurs  within  this  diffusion  mantle  because  of  the  extreme 
temperature  gradient.  Exothermic  decomposition  may  change 
tho  mantle  temperature  gradient  to  such  an  extent  that  a 
decomposition  flame  appears  within  the  boundary  layer. 

This  flame  has  the  effect  of  maintaining  high  vaporization 
rates  even  under  low  convection  (flow)  conditions.  Thus 
the  physical  two-flame  model  presented  here  determines 
higher  combustion  rates  near  the  injector  face  than  does 
the  single-flame  model  of  Reference  1  and  10. 
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2. 


The  Model 


Experiments  and  a  literature  search  indicated  that 
decomposition  and  oxidation  flame  fronts  coexisted  in 
the  combustion  of  N2H4  type  fuels.  Thus  a  model  which  con¬ 
siders  two  different  burning  mechanisms,  depending  upon 
droplet  Reynolds  number,  was  developed.  Such  a  model  is 
compatible  with  the  single-mechanism  system  of  equations 
developed  in  Reference  1. 

a.  Oxidation  Model  (Single-Flame  Regime). 

This  model  incorporates  features  of  two  other  models 
in  the  literature  (Refs.  1,  10).  The  treatment  of  vapori¬ 
zation  which  allows  for  heating  of  the  drops  is  essentially 
that  of  Priem  and  Heidmann  (Ref.  1)  while  the  treatments  of 
drop  ballistics  and  chamber  gas  dynamics  are  essentially 
those  developed  b/  Rocketdyne  (Ref.  10).  The  combustion- flow 
process  is  considered  to  be  one-dimensional.  The  model  re¬ 
presents  the  continuous  distribution  of  drop  sizes  by  an 
arbitrary  number  of  size  groups  as  do  both  references.  A 
logari thmi conormal  distribution  has  been  used  to  describe 
the  spray  pattern  of  the  motor  of  interest;  however,  any  * 

type  distribution  may  be  substituted  into  the  droplet  size 
subroutine. 

(1)  System  Equations.  The  equations  describing 
combustion,  assuming  vaporization  is  rate  controlling, 
are  presented  in  this  section.  These  equations  are  mani¬ 
pulated  so  that  simultaneous  computer  solution  is  possible. 

The  existence  of  a  decomposition  front,  within  the  boundary 
layer  mantle,  is  consistent  with  these  equations  except  that 
a  different  mass  vaporization  rate  must  be  defined  as  dis¬ 
cussed  in  Section  Ill-2,b.(2). 

(a)  Droplet  Ballistics.  A  Force-momentum 
balance  on  the  drop  yields  the  equation 

dv  3Ca°  8  I u-vl  (u-v) 

dx  "  8Pt  r  v  (II  I- 1) * 

*  Equations  following  by  an  asterisk  are  applied  individually 
to  all  drop  sizes  of  both  oxidizer  and  fuel. 
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where  is  defined  as  in  reference  10. 


C.  ■  27  Re 

d  g 


•  0.84 


«  0.271  Re 


0.217 
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0<  Re  <80 
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80<Re  <10 
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(III-2)* 


Reg  “ 


2rA  vp, 
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104<Re. 


( 1 1 1  -  3)  * 


(b)  Droplet  Mass.  The  rata  of  evaporation 
from  a  spherical  drop  of  propellant  is  given  by  the  equation 
(Ref.  1). 

2ir  D„r  «Nu 


'a*  a  "“a 


RT 


(III  —  4)* 


where 


ln  l£rp“J»  and 


( 1 1 1  -  5)  • 


the  Nusselt  nuaber  for  mass  transfer  is  determined  from 
Ranz  and  Marshall's  correlation  (Ref.  11). 


Nu_  ■  2  ♦  0.6  Re  1/2  Sc1/3 

m  m 


(  1 1 1  -  6  )  * 


Thus  the  change  in  droplet  mass  with  distance  is  given  by 

6 


da  a  da  dt  a  w 

17  17  17  v  (1 1 1  -  7 )  * 

(c)  Heat  Transfer  to  Droplet.  A  heat  balance  on 
a  single  droplet  yields  the  result  that  the  rate  of  accumula¬ 
tion  of  en.rgy  must  equal  the  net  energy  transferred  to  the 
drop  at  any  tiae. 


•Equations  followed  by  an  asterisk  are  applied  individually 
to  all  drop  sizes  of  botl.  oxidizer  and  fuel. 
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Thus  , 


«  c„  ill  -  qv  -  .  w,l  •  Qr  (III-S) 

”  1  Q  fl  2 

The  term  QR  represents  the  radiant  heat  transfer  to  the  drop 
and  has  been  omitted  from  this  program;  although,  it  could 
be  included  if  values  for  emissivity  justify  it. 

The  third  term  on  the  right  hand  side  of  (III-8)  represents 
the  kinetic  energy  imparted  to  the  vapor  leaving  the  surface. 

This  term  is  generally  negligible  but  was  included  in  the 
program  since  it  can  be  important  for  small  radii  and/or  for 
drops  near  their  critical  temperature  (for  which  X  approaches 
zero) . 

The  term  represents  the  heat  transferred  to  the  vaporizing 

drop  by  convection  with  mass  transfer.  This  is  defined  in 
terms  of  a  heat  transfer  coefficient  as 


qv  -  4*r2h  (Tg-T4)  (HI-9) 

Where  the  heat  transfer  coefficient  with  mass  transfer,  h, 
is  related  to  the  coefficient  without  mass  transfer  by 


h 


and  z  ■  WCP« »  wCP».p 

4*r2h  2»r  Nu^ 

i 

Again  the  Nusselt  number  for  heat  transfer  has  been  defined 
in  the  work  of  Ranz  and  Marshall  (Ref.  11)  to  be 

Nuh  .  2  ♦  0.6  Re^2  Pr1/3 
Substituting  into  (III-9) 

q  -  "  CP,,m  <VT*> 

(«l  -1) 


(III-10) 

(III-ll)* 


(III  —  12) * 


(III-13)* 


Finally  the  heat  balance  ( 1 1 1 - 8 )  yields 

art  _  qv-«x  - 

dx  "  - — - 

m  Cp  £v 


(II 1-14) * 
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(d)  Chamber  Pressure.  A  force  balance  on  an 
incremental  chamber  length  yields 


d(PAc)  dI 

“3*  "’77 


where  I  ■  P g  u  ^c  ♦In 

g 

The  continuity  equation  is  expressed  as 

Pg  u  Ac  *  En  (mj  -■) 


a  v 


(III  - 15) 


(III-16) 


(111-17) 


Use  of  equations  (JII-16)  and(III-17)  in  equation  (III- 15)  yields 

£  ■  tIn(”J'm)  £  * En  [■  £ '  tu-v)  £]• p  I  (m-18) 


Thus  the  pressure  gradient  has  been  expressed  in  terms  of 
velocity  gradient.  However,  another  independent  linear 
relationship  between  du  and  dP  must  be  developed  so  that  dP 

77  17  '  .  77 

may  be  computed  without  iteration. 

The  following  equations  are  necessary  for  this  development. 


♦ 

>  ^min  j’B^o 

(III  —  19) 

[En(minj-m) ]f 

T 

g 

■  T  [l-  u2(Y~1)H«l 

(III-20) 

L  2y  RT0  g  J 

«  P  Mg 

(III-21) 

The  analysis  assumes  the  ability  to  evaluate  properties  of  the 
combustion  product  gas  (Tg0,  M  ,  y  ,  k  ,  and  v_)  as  functions 
of  chamber  pressure  and  mixture  ratio  and  propellant  properties 
(Pa*  ka.  pt  »  y *  cDf»  »nd  Cpa)  as  functions  of 


tern; erature . 


a*  *1  •  T*  ''pi* 

Furthermore,  the  diffusion  coefficients,  D,  will  have  to  be 
determined  functions  of  mixture  ratio,  pressure,  and  tempera 
tare. 
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After  considerab  c  algebra  involving  equations  (II  1-17,  19, 

20,  and  21)  and  their  derivatives,  an  independent  linear 
relationship  between  du  and  dP  can  be  developed  in  the  form 

dx  Tx 

37  ‘  F  *  G  37  (In-22) 

where  F  and  G  are  functions  of  quantities  which  are  known 
at  any  step  in  the  numerical  integration  (e.g.,  m,v,dm,  dv 

17  Tx 

P,  u,---),  of  the  combustion  product  gas  properties  and  of 
the  derivatives  of  those  properties  with  respect  to  local 
mixture  ratio,  ♦  . 


If  equation  (111-18)  is  considered  to  have  the  form 


■  h  ♦  j  in 

dx  dx 

equations  (III-22  and  23)  result  in 


( 1 1 1  -  2  3) 


dP  -  11  4  JF 
I  -  JG 


(III  —  24) 


Thus,  equation  ^xII-24)  allows  dP  to  be  computed  directly 
in  this  program.  Tx 


(2)  Physical  Properties  in  the  Vaporization  Mantle. 


(a)  The  NTO  System,  The  nitrogen  tetroxidc  system 
exhibits  complex  changes  with  very  short  equilibrium  times 
(order  of  microseconds).  The  system  exists  in  equilibrium 
in  the  following  form 


N2°4  *■  2N02  ♦  2N0  +  °2 


The  extent  of  either  of  the  two  decomposition  reactions  is 
a  function  of  the  system  temperature  and  pressure.  The 
actual  function  is  not  well  defined  by  thermodynamic  data 
because  N2O4  and  NO2  exist  in  equilibrium  at  the  mixture 
critical  point.  Thus  it  is  impossible  to  separately  determine 
either's  critical  properties. 

For  the  temperature  anu  pressure  ranges  to  be  considered  in 
this  work  it  was  possible  to  postulate  an  approximate  physical 
model.  Utilizing  such  a  model,  liquid  and  vapor  composition 
and  thus  such  properties  as  thermal  conductivity,  heat  capa¬ 
city,  and  heat  of  vaporization  could  then  be  computed. 

As  shown  in  Figure  9,  the  entering  liquid  droplet  is  assumed 
t  )e  N2O4.  However,  as  the  N2O4  evaporates  it  goes  to 
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The  assumption  making  possible  the  inclusion  of 
heat  of  reaction  with  heat  of  vaporization  is 

6r  <  10%  of  (R  -  r  ). 


FIGURE  9.  Nitrogan-Tatroxide  Droplat  Vaporization  Modal. 
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approximately  50%  NO,.  Then  since  the  temperature  gradient  is 
very  steep,  (Figure  9)  the  vapor  will  be  100%  NO2  before  it 
crosses  even  10%  of  the  diffusion  mantle.  The  approximation 
that  NO2  is  the  vaporizing  species  was  thus  made.  Using 
this  approximation  it  was  possible  to  account  for  chemical 
change  by  adding  the  heat  of  decomposition  to  the  heat  of 
vaporization  for  N2O4.  The  vapor  was  assumed  to  consist  of 
NO2  for  the  purpose  of  physical  property  calculation.  The 
second  decomposition  2NO2  2  2NO  ♦  O2  was  not  considered  in 
this  approximate  model. 

(b)  The  MM1I  System.  Monomethy  lhydrazine  can 
support  a  decomposition  front  in  the  vapor  mantle  before 
undergoing  oxidation.  The  occurrence  of  such  phenomena  not 
only  changes  such  properties  as  thermal  conductivity,  heat 
capacity,  density  and  diffusivity  but  also  necessitates 
postulation  of  a  decomposition  model.  A  physically  realistic 
model  for  va.orization  with  decomposition  is  discussed  in 
the  next  section. 

b.  Decomposition  Model  (Two-flame  Regime). 

(1)  Introduction.  The  equations  presented  in 
Section  2. a.,  describe  a  complex  two-phase  gasdynamics 
system.  The  vaporization  description  is  one  part  of  this 
coupled  system.  This  part  was  modified,  so  that  steady- 
state  combustion  with  vtror  phase  decomposition  of  a 
thermally  unstable  fuel  could  be  described.  The  modified 
combustion  package  was  then  incorporated  into  the  system 
description. 

Dynamic  Science  Corporation  (Ref.  12)  has  completed  experi¬ 
ments  which  indicate  that  a  single  diffusion  flame  front 
model  is  not  realistic  for  such  monopropellant  systems  as 
hydrazine  and  monomethy lhydrazine.  As  shown  in  Figures  10 
and  11,  two  flame  fronts  appear  in  the  combustion  of  hydrazine 
and  MMH  in  nitrogen  tetroxide.  The  presence  of  these  two 
flame  fronts  has  been  observed,  but  the  exact  nature  of 
these  fronts  has  not  been  fully  determined. 

This  section  reviews  experimental  and  analytical  studies 
of  the  MMH  type  system.  A  realistic  two  flame  model  is 
postulated  along  with  its  required  assumptions.  Equations 
are  derived  for  this  two  flame  decomposition-oxidation  model 
and  the  method  of  solution  is  explained.  The  results  of  this 
work  are: 


i)Speci fication  of  a  realistic  steady-state 
vaporization  model  which  includes  vapor 
phase  decomposition, 
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Note:  The  liquid 
drop  is  within  the 
wire  support. 


Hydrazine  burning  in  still  nitrogen 
tetroxide,  at  one  atmosphere 


Photograph  of  Expariaantal  Burning  Drop 
Showing  Two  Flaae  Fronts  (Raferenca  12) 
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FIGURE  10. 


Note:  The  liquid 
drop  is  within  the 
wire  support. 


MMH  Burning  in  still  nitrogen 
tetroxide,  at  one  atmosphere 


I 


FIGURE  11.  Photograph  of  Experimental  Burning  Drop 
Showing  Two  Flame  Fronts  (Reference  12) 
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ii)  Development  of  the  technique  used  to 
solve  these  equations. 

iii)  Calculation  of  the  spatial  combustion  pro¬ 
perties  with  vapor  phase  decomposition. 

The  vaporization-diffusion  model  with  a  single  stoichiometric 
diffusion  flame  front  is  shown  in  Figure  12b.  For  bipropellant 
combustion  the  outside  boundary  of  the  vapor  mantle  determines 
the  position  of  the  stoichiometric  diffusion  flame  front. 
However,  monopropellants  may  support  a  decomposition  front 
within  this  diffusion  mantle.  The  decomposition  front  is 
similar  to  a  premixed  gas  flame  front  and  cannot  be  analyzed 
as  a  diffusion  flame.  When  a  decomposition  front  occurs 
within  the  fuel  diffusion  mantle  of  an  oxidizer-rich  propel¬ 
lant  combustor,  two  flame  fronts  will  appear.  The  inside 
one  is  a  decomposition  front  while  the  outer  flame  is  an 
oxidizer  diffusion  flame  front. 

The  postulated  physical  model  involves  two  regions: 

Region  I,  (Figure  12a)  where  the  two-flame  decomposition 
mechanism  is  rate  controlling,  and 

Region  II,  (Figure  12b)  where  the  single-flame,  laminar 
boundary- layer  oxidation  mechanism  is  rate  controlling. 

The  equations  applicable  to  Regions  I  and  II  were  used  in 
specific  calculations  for  the  MMH-NTO  system.  The  calcu¬ 
lated  burning  rates  for  a  distribution  of  MMH  drops  were 
determined  and  used  to  calculate  the  fraction  combusted 
and  the  C*  efficiency.  The  C*  efficiencies  for  several 
engines  were  finally  compared  with  experimental  values 
determined  by  Priem  (Ref.  1). 

(2)  Equations  and  Calculations.  The  existence  of 
either  region  I  or  II  depends  upon  the  relative  velocity 
between  the  liquid  drop  and  the  product  gas.  Region  II 
corresponds  to  a  large  relative  velocity  (and  small  boundary 
layer),  while  region  I  corresponds  to  a  small  relative 
velocity  (and  relatively  stagnate  conditions).  The  applica¬ 
bility  of  the  equations  of  either  region  are  thus  dependent 
upon  the  criterion  that  either  the  decomposition  thickness  is 
within  the  boundary  layer,  rt<B  (region  I),  or  the  boundary 
layer  is  less  than  the  decomposition  thickness  rt>B  (region  II). 
The  equations,  which  allow  calculation  of  r*  and  B,  were 
developed  and  used  simultaneously  to  determine  applicable 
regions . 
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Figure  12a.  Decomposition  Front  Within  Diffusion  Mantle,  Regime  I 


upon  Nusselt  number  correlation 


Finite  decomposition 
gaseous  front  where 
r.  is  independent  of 
Flow 


Diffusion  flame  front 
(oxidation) 


Figure  12b.  Laminar  Diffusion  Model  when  rt5»  B,  Regime  n 


Decomposition  and 
diffusion  flame  front 
together  when  rt»  B 

Note:  B  decreases 
with  velocity  while  rt 
is  relatively  independent 
of  velocity 


Figure  jj.  Two-Regime  Model  of  Decomposition -Front  Position, 

Adiabatic  Oxidation  Temperature  the  same  in  both  cases. 
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Decomposition  Equations  (Region  I) 


The  theoretical  decomposition-flame  model  of  Tarifa  and 
Notario  (Ref.  13)  was  extended  to  include  the  influence  of 
adiabatic  oxidation  flame  temperatures.  The  equations 
applicable  for  small  drops  in  this  region  are  not  derived 
here;  but  are  presented  from  the  literature  (Ref.  13  )  so  that 
different  liquid  properties  may  be  substituted  in  the  program. 


The  decomposition  film  thickness  is  given  in  reference  13  as, 


X. 

-1  r , 
X* 


(III-2S) 


And  the  decomposition  evaporation  rate  constant  (in*/§ec), 

[defined  by  K  ■  fLil—L]  given  in  reference  13  bv 

dt 


Y  m 

2k  X 

(III-26) 

N  ■ 

Cp  Pl 

where 

xs  “ 

lr. 

[®m  ♦  X* 

1  et  -C 

(1 1 1 -27) 

X* 

3 

"  4/r 

rln(e./9j)  i  ■  a. 

1(8.  -  exp  2(e.  -  \) 

r  (II 1-28) 

and  0 

■ 

<lr 

(T  -  T.+  i_  ) 

4  c 

cp 

(I I I -29) 

e 

•  0 

m 

for  T  ■  Tm  (adiabatic  flame  temperature) 

9 

•  ei 

for  T  ■  T^  (liquid  temperature) 

e 

"  eo 

for  T  ■  0 

•a 

-  C_E 

lq7 

(III -30) 

\ 

•  B 

h  bt  TiCp\n 

(III-31) 

k  V  qr  i 

where 

(Ref. 

B  and 
13). 

n  are 

constants  of  the  reaction  velocity 

equation 

In  order  to  use  the  model  proposed,  the  decomposition  thickness 
and  the  evaporation  rate  constant  must  be  calculated  and  used 
as  input  to  the  program.  The  decomposition-flame  film  thickness 
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and  the  evaporation  rate  constant  were  determined  for 
several  adiabatic  oxidation  temperatures  and  arc  shown  in 
Figures  13and  14for  the  MM1I-NT0  system.  An  adiabatic 
decomposition  flame  temperature  for  a  monoprope 1 lant 
corresponds  to  Qm  ■  1.0;  while  =  1.54  corresponds  to  an 
adiabatic  oxidation  flame  temperature  for  the  MMH-N'TO  bi- 
propellant  system.  The  calculated  curves  of  Figuresu  and 
are  thus  extensions  of  the  work  of  reference  13  to  two 
flame  front  conditions. 


(3)  Calculation  Procedure  with  Two-Flame  Model. 

The  decomposition  thickness  and  evaporation  rate  curves 
of  Figures  13andi4and  burning  rate  curves  similar,  to 
Figure  15  were  used  simultaneously  to  compute  burning 
rates  for  the  two  regions.  The  values  for  the  single 
flame  model  are  shown  in  dashed  lines  (Fig. 15)  far  comparison, 

Burning  rate  curves  similar  to  Figure  15  are  computed  by 
the  following  steps  for  each  group. 

t 

(a)  Determine  regions  by  the  criteria. 


rt  <  B 
rt  >  B 


region  I 
region  II 


Where  B  is  the  boundary  layer  thickness  shown  in  Figure  15 
and  rt  is  the  decomposition  film  thickness  shown  in 
Figure  13  (B  and  rt  are  input  to  the  program). 

(b)  By  this  criteria  the  burning  rate  of 
Figure  15was  initially  determined  by  the  mechanism  of 
region  II. 


of  region  I. 


(c)  At  rt  ■  B,  the  mechanism  shifted  to  that 


(d)  The  burning  rate  for  region  II  was  now  deter¬ 
mined  by  the  K  shown  in  Figure  14  and  the  rate  equation: 


3  rK 

7  “ 

rov 
(e) 

that  of  region  II . 


fraction 

inch 


(III-32) 


At  rt  ■  B  the  mechanism  shifted  back  to 


(f)  If  the  drop  had  not  yet  reached  its  wet  bulb 
temperature,  while  B  >r^,  the  distance  to  do  so  was  calculated 
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Decomposition  Film  Thickness,  rt ,  (microns) 


125 


0  100  200  300 
Drop  radius,  r,  (microns) 


FIGURE  13.  Film  Thickness  For  Decomposition  Flame 
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Evaporation  Constant 


Film  thickness,  B  (microns)  Combustion  rate,  infraction/ inch) 


l 


FIGURE  15 .  Combustion  Rate,  Initial  Drop  Size  49.  4v 
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by  a  heat  balance  which  resulted  in  the  equation 


Ax  °  P  w .  b , 

2*  <\r  pr  k 

(g)  The  resulting  burning  rate  curves,  shown 
in  Figure  IS  arc  continually  integrated  by  the  program  to 
determine  the  cumulative  fraction  of  fuel  vaporized.  This 
curve  is  shown  in  Figure  16. 

(4)  Characteristic  Results  Using  Two-Flame  Model. 

Burning  Rat es .  The  burning  rates  computed  for 
the  decomposition  model  are  substantially  higher  than  those 
computed  using  the  single  flame  front  model. 

The  maximum  burning  rates  under  either  the  single  or  double 
flame  models  are  essentially  the  same  (once  wet  bulb  tempera¬ 
ture  is  reached)  because  adiabatic  conditions  are  assumed  and 
the  same  amount  of  total  enthalpy  release  is  involved  in 
both  mechanisms. 


The  zero  relative  velocity  point  has  a  less  pronounced  effect 
upon  the  burning  rate  because  the  decomposition  flame  now  de¬ 
termines  the  minimum  rate  at  stagnation  conditions. 

Fraction  Vaporized.  The  fraction  vaporized 
under  the  two  flame  model  is  higher  than  under  the  single 
flame  model  as  shown  in  Figure  16,  At  one  inch  chamber  length 
the  fraction  is  approximately  100%  greater  for  the  two-fl&mc 
model  than  for  the  single-flame  model. 


C*  Efficiency ,  nc»t  j^e  experimental  data 
shown  in  Figure  17were  used  by  Priem  and  He i draann ( Ref .  1) 
to  compare  experimental  and  calculated  hydrazine  efficiencies. 
This  comparison  for  a  large  number  of  propellant  combinations 
is  shown  in  Figure  18, 


The  hydrazine  fraction  vaporized,^  ,  (Figurei3)  is  vaporiza¬ 
tion  rate  controlling  with  the  gaseous  oxygen  data  of  Fig.  l?.. 
Thus,  the  computed  was  used  with  6  ■  1.0  for  gaseous  oxygen 
to  compute  the  efficiency  for  hydrazine  -  COX  according  to  the 
formula 


C* 


(C»)  O/Jf 
(C*)  O/F 


♦  1 


(III-33) 


3S 


t 
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C*  EFFICIENCY 


.8 


Chamber  Length,  Inches 


IGURE  17,  Experimental  Performance  of  Hydrazine 
Propollanr,  Combinations  (  Ref.  1). 

37 


FIGURE  18.  Coaparison  of  Expariaental  and  Calculated  Efficiencies  (Ref.  1) 
Showing  Hydrazine  Efficiencies  computed  by  the  two-flaae  aodel 


The  efficiencies  computed,  using  the  two-flame  model,  are 
in  good  agreement  with  the  experimental  efficiencies  as 
shown  in  Figure  is. 

The  two-flame,  vapor  decomposition  model  of  droplet  combus¬ 
tion  in  hydrazine  type  fuel  systems  allows  accurate  pre¬ 
diction  of  the  droplet  history  in  a  combustion  chamber. 

This  information  is  necessary  for  both  efficiency  and 
instability  calculations.  This  model  therefore,  represents 
a  significant  advance  in  the  state-of-the-art  of  analytical 
chamber  design. 

3.  The  Stability  Drop  Size 

The  model  used  for  the  analysis  of  stability  treats  drop¬ 
lets  as  if  they  were  all  of  a  single  size.  It  is  desirable  to 
select  a  drop  size  which  represents  the  vaporization  properties 
of  the  spray  distribution  as  accurately  as  possible. 

The  vaporizTTTon  rate  of  a  spray  made  up  of  a  finite  number 
of  drop  size  groups  (as  used  in  the  steady-state  combustion 
model)  is  given  by 


u  ■ 

I  ci®iri  (2  ♦  hr^ 

1/2  h  V 

(111-34)- 

where 

Bi  ’ 

p 

p 

(  1 1 1 -  35) 

p-p.i 

♦i  ' 

c 

• 

< 

K> 

( 1 1 1  -  36) 

Ci  ' 

(W‘«i>,/2 

(III-37) 

K 

.6  /7  scl/5 

(III-38) 

Combustion  instability  is  related  to  changes  in  u 
by  disturbances  in  the  state  of  the  gas  (p,T)  and 
relative  velocity  term,  4  .  If  we  preserve  w,)w, 

from  a  multigroup  representation  to  a  single-group 
following  equations  are  obtained: 

caused 
in  the 

3w,and  du 

7p  TT 

model,  the 

39 
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w 


Equation 

c*r»B*(2*kril/24tCt)' 
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EciriBi*2*kri  M*) 


(111-39) 
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3/2  3/2 

csrsBs^s  *  J  ciri®iCi 


(in-40) 


is  c»ril2  ff  ♦  k  i-p  •  «.  “  )) 

-  J'ini2  ff  *  k  n/2n<q  if  *\  If  )l  (in-4i) 

9w  c  r  f2  ^  4  k  f  r  ££  ♦  a  )  1 

IT  csV2  If  *  Bs  It  n 

-  2  «i»iH  |f  ♦  kn/2»i(«i  If  If  >1  c,!I-42> 


In  equations  (111-41  and  42)  it  has  been  assumed  that  the  spray 
of  cnly  one  propellant  is  controlling  the  stability  and  that 
0  and  £  can  be  represented  by  functions  which  are  linear  in 
p  and  T  so  that  the  partial  derivatives  are  constants  determined 
from  the  propellant  properties. 


Partial  solutions  to  equations  (III-39  -  42)  can  be  obtained 
under  certain  assumptions.  If  we  assume  that  ail  +  are  large, 
thi.t  96 

Ip 

obtain 


■  30 

If 

■  0  and,  as  a  result, 

that  all  Bj  are  equal,  we 

3/2 

a  a 

Eci  ri  Mi 

?s 

J/2. 

(III-43) 

I  Vi  'i 

Cs  -  £ciri3/2Ci4i 

- m — 

Eciri  U 


( 1 1 1  -  44 ) 


The  assumptions  made  to  arrive  at  equations  (III-43  and  44) 
are  equivalent  to  those  w*ue  by  Priem  in  his  stability  model 
(Ref.  4),  namely  that  changes  in  vaporization  rate  are 
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attributed  only  to  changes  in  the  Nusselt  number  and  that 
the  Nusselt  number  is  much  larger  than  2. 

On  the  other  hand,  if  interest  is  focused  on  regions  of  small 
velocity  difference  where  stability  criteria  seem  to  be  most 
critical,  we  seek  a  solution  for  small  The  results  are 


8s 


rciriBi 

ICiTi 
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1/2 

s 


*s 


£c  r  3/2  8  £ 
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l  cir.6i 


(111-45) 

(III -46) 


The  quantity  of  pm/i'm  is  much  less  sensitive  to  drop  temperature 
than  8.  If  the  are  considered  to  be  nearly  equal,  equation 
(III-46)  can  reduce  to 


3/2 

Eciri  Bj 
Eciri®i 


(II 1-47) 


In  the  region  between  the  injector  and  the  point  where  drop 
velocity  and  gas  velocity  are  equal,  small  drops  are  especially 
important  for  two  reasons;  1)  The  absolute  velocities  of 
small  drops  are  lower  than  those  of  larger  drops,  resulting 
in  increased  concentration,  c,  of  small  drops.  2)  The  small 
drops  heat  up  faster  and  attain  high  vaporisation  rates  due 
to  high  value  of  8  which  are  quite  sensitive  to  the  large 
increases  in  vapor  pressure  which  accompany  the  rise  in  drop 
temperature . 

This  analysis  and  results  from  the  steady-state  combustion 
program  indicate  that  combustion  instability  may  be  governed 
by  the  properties  of  the  smaller  drops  in  a  spray  rather 
than  by  those  of  the  drop  of  mass  mean  radius. 
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IV.  COMBUSTION  INSTABILITY  MODEL 


1 .  Introduction. 

While  there  has  been  a  great  deal  of  research  devoted 
to  create  linear  models  of  the  instability  phenomena,  little 
has  been  done  on  nonlinear  models,  in  spite  of  the  fact  that 
combustion  instability  in  liquid  rocket  engines  is  z  nonlinear 
phenomena  (Ref.  8  and  9) .  The  fundamental  conservation  equations 
i.e.,  mass,  momentum,  energy,  and  state  will  result  in  a  non¬ 
linear  model,  by  their  very  nature.  Because  of  the  complexity 
of  the  solution  of  these  nonlinear  partial  differential  equa¬ 
tions,  various  authors,  among  them,  Crocco  (Ref.  14)  and 
Priem  (Ref.  1),  have  attempted  to  simplify  them  before  solu¬ 
tion  with  various  assumptions.  The  conservation  equations 
in  their  most  general  form  consist  of  four  independent  variables 
(6,  r,  z,  t) ,  making  a  closed-form  solution  impractical.  With 
the  advent  of  the  high-speed  digital  cjmputer  and  advanced 
numerical  techniques,  solving  these  equations  becomes  possible. 
However,  solution  with  four  independent  variables  is  a  para¬ 
mount  task  and  would  exceed  the  practical  cost  of  computer 
runs.  Priem  has  solved  these  equations  with  two  independent 
variables  (0,t)  with  simplifying  assumptions.  In  this  approach, 
functional  dependencies  are  determined  in  every  case  from  the 
basic  equations  for  the  phenomena,  and  simplifying  assumptions 
are  introduced  that  are  considered  reasonable  in  light  of  the 
combustion  processes  in  a  liquid  rocket  engine.  The  models 
advanced  by  Dynamic  Science  Corporation,  Crocco  and  Priem 
utilize  the  same  theory  derived  from  the  fundamental  conser¬ 
vation  equations,  but  with  different  assumptions. 

Dynamic  Science  Corporation  feels  that  combustion  insta¬ 
bility  in  liquid  rocket  engines  is  a  nonlinear  phenomena,  in 
that*  1)  stability  depends  on  the  disturbance;  2)  there  is  a 
limiting  amplitude  of  oscillations;  and  3)  the  wave  forms 
become  steep  fronted  or  nonsinuisoidal.  In  addition,  the  non¬ 
linear  model  enables  the  combustion  process  to  be  described 
from  the  basic  equations  of  the  phenomena  without  the  intro¬ 
duction  of  "empirical"  or  "intuitive"  constants.  Linear  models 
cannot  relate  engine  parameters  to  threshold  disturbance. 

Either  a  linear  system  is  stable  or  unstable  regardless  of  the 
disturbance.  Such  linear  models  are  derived  from  the  nonlinear 
conservation  equations  by  assuming  that  all  variations  are  small. 
This  assumption  is  necessary  if  a  closed  form  solution  is  de¬ 
sired  since  mathematical  techniques  are  not  well  developed  for 
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solving  simultaneous  sets  of  nonlinear  partial  differential 
equations . 

Dynamic  Science  Corporation  model  is  based  on  the  numeri¬ 
cal  solution  of  the  conservation  equations  with  mass-addition. 

The  mass- addition  or  combustion  rate  is  described  by  a  quasi¬ 
steady-state  expression  for  vaporization.  This  assumes  that 
all  the  combustion  processes,  i.e.,  atomization,  mixing,  and 
kinetics  are  rapid  compared  to  the  vaporization  process.  The 
quasi-steady-state  assumption  implies  that  the  droplet  spray 
instantly  responds  to  a  change  in  the  local  environmental 
conditions . 

2 .  Theory . 

An  annular  section  of  a  combustion  chamber  is  chosen  to 
predict  the  response  of  the  spray  to  a  disturbance.  The  use 
of  an  annular  section  is  dictated  by  the  limitations  of  present 
computers  as  to  storage  and  speed  in  producing  useful  results. 

An  annular  section  enables  the  variations  of  the  dependent 
variables  in  the  tangential  direction.  Such  a  model  enables 
the  description  of  a  tangential  disturbance  propagated  through¬ 
out  the  spray  field.  As  shown  in  Figure  1  the  annular  section 
has  a  very  small  length  Az  and  thickness  Ar.  Physically  this 
means  that  there  are  no  gradients  in  the  r  and  z  direction 
in  the  annulus.  It  is  assumed  that  the  steady-state  combustion 
is  uniform  around  the  annulus.  However,  the  transient  combus¬ 
tion  rate  is  a  function  of  the  local  environment  and  varies 
around  the  annulus.  The  assumption  that  steady-state  properties 
are  not  a  function  of  6,  implies  that  the  propellants  are  in¬ 
jected  uniformly, at  the  injector  face, at  that  particular  radius 
of  injection.  The  annular  model  can  be  extended  to  injection 
profiles  which  are  distributed  across  the  injector  face  by  use 
of  the  "nodal  method"  and  a  quasi  two-dimensional  steady-state 
model.  While  the  annular  model  is  one-dimensional,  in  order  to  de¬ 
scribe  the  convective  flow  out  of  the  annulus,  the  z  derivatives 
are  treated  in  a  gross  way  by  assuming  that  the  dependent 
variables  are  constant  through  the  annulus  but  the  derivatives 
vary  with  time.  The  model  does  not  treat  the  z  direction  in  a 
purely  two-dimensional  manner  since  the  z  derivatives  are 
lumped  in  the  6  direction  by  an  integration  assuming  mass, 
momentum,  and  energy  is  constant  in  the  annulus  as  a  function 
of  time.  The  derivation  of  the  transport  or  field  equations 
is  similar  to  the  method  of  Priem-Guentert  (Ref.  1)  and  will 
be  repeated  here  for  completeness  of  the  nonlinear  model. 

The  basic  transport  equations  for  two-phase  flow  have 
been  developed  by  Bird,  Stewart,  and  Lightfoot  (Ref.  15): 
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Equation  of  Continuity 
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After  assuming  the  liquid  velocity  and  temperature  to  be 
constant  in  the  annulus  and  further  manipulation  of  the  equa¬ 
tions  IV- 1 ,  2,  3,  we  obtain  the  transport  equations  for  two 
phase  flow  with  mass  addition. 

The  conservation  equations  appear  in  the  following  form: 
Continuity 
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The  above  equations  can  be  nondimcnsionali  zed  by  means  of 
the  following  transforms: 
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Apply  the  transforms  to  the  continuity  equation 
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After  further  manipulation  the  resulting  equations  are: 
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Energy 


I 


By  assuming  that  the  gas  is  perfect  and  flow  is  one-dimensional 
and  isentropic: 

»o2  *  8  tfrTo  ( I V- 2  1 ) 


then  equations  IV-18,  19,  and  20  can  be  reducted  to 
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Equations  (IV-33)  is  similar  to  Damkohler's  group  bassd 
on  the  speed  of  sound.  It  can  be  thought  of  as  the  ratio  of 
the  wave  time  around  the  annulus  at  the  speed  of  sound  to  the 
combustion  time.  L  is  defined  by  Priem  as  the  burning  rate 
parameter.  Equation  (IV-35)  is  similar  to  the  reciprocal  of 
a  Reynolds  number  based  on  the  speed  of  sound.  J  is  defined  as 
the  viscous  dissipation  parameter  and  is  a  measure  of  the 
energy  lost  by  viscous  forces.  By  substitution  of  Equations 
(IV-33)  and  (IV-3S)  the  following  transport  equations  are 
obtained : 
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(  I V-  36  ) 
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Using  the  appropriate  cylindrical  coordinate  transforms 
outlined  in  Bird,  Stewart,  and  Lightfoot,  (Ref.  15),  equations 
(IV-37-39)  can  be  simplified  to  describe  an  annular  chamber 
section.  If  it  is  assumed  that  all  dependent  variables  and 
derivatives  of  dependent  variables  in  the  radial  direction 
arc  zero,  that  all  axial  dependent  variables  do  not  vary  with 
angular  position,  and  that  all  second  derivatives  in  the  axial 
direction  arc  zero,  the  following  transport  equations  result: 
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Momentum  (in  the  axial  or  z-dlrection) : 
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To  obtain  the  lumped  derivative*  of  the  axial  dependent 
variable  to  account  for  convective  flow  through  the  annulus, 
it  is  assumed  that  conservation  of  mass,  momentum,  and  energy 
must  apply  for  the  entire  annular  system  as  well  as  the 
incremental  volumes  described  by  equations  (1V-41-43). 
Following  the  Priem-Guentert  (Ref.  4)  assumptions  to  determine 
the  derivatives  in  the  axial  or  z  -  direction  the  conservation 
equations  are  integrated  through  the  entire  volume  of  the 
annulus.  Since  it  is  assumed  that  none  of  the  terms  vary  in 
the  axial  or  radial  direction,  the  equations  only  have  to  be 
integrated  in  the  e-direction.  Further,  if  it  is  assumed 
that  the  total  mass,  momentum,  and  energy  in  the  annulus  does 
not  vary  with  time.  It  is  also  assumed  that  all  dependent 
variables  are  zero  in  the  radial  direction  and  constant  across 
Az.  All  axial  derivatives  do  not  vary  with  angular  position 
and  are  thus  lumped  or  averaged  by  the  integration.  Since 
the  annulus  is  closed  in  the  angular  direction  there  can  be 
no  net  loss  or  gain  of  mass,  momentum,  and  e  ergy  in  that 
direction.  These  assumptions  result  in  the  following  integral 
equations  for  the  axial  derivatives: 
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To  define  the  annular  model  an  expression  for  the  mass 
addition  is  required.  In  a  liquid  rocket  engine  four  signi¬ 
ficant  processes  determine  the  mass  addition  or  combustion 
rate:  1)  atomization;  2)  vaporization;  3)  mixing;  and  4) 
kinetics.  It  has  been  shown  by  many  authors  that  for  a  well 
mixed  chamber  at  typical  pressures  and  temperatures  encountered 
in  a  liquid  rocket  engine  that  atomization,  mixing,  and  kinetics 
are  much  faster  processes  than  vaporization.  Thus  it  may  be 
concluded  that  the  combustion  rate  can  be  assumed  to  be  con¬ 
trolled  by  vaporization.  Further  in  a  well  mixed  and  atomized 
spray  the  combustion  rate  is  controlled  by  droplet  vaporization. 
Since  high  performance  liquid  rocket  engine  systems  are 
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controlled  by  droplet  vaporization,  the  instability  model 
describes  the  response  of  a  vaporizing  droplet  spray. 

With  the  above  description  of  the  droplet  combustion  model 
shown  in  Figure  19,  a  vaporization  rate  equation  may  be 
derived. 

The  vaporization  rate  equation  may  be  derived  analyti¬ 
cally  (Ref.  1)  in  terms 


However,  the  state-of-the-art  does  not  allow  for  analy¬ 
tical  calculation  of  the  flow  dependent  film  thickness,  B . 

Thus,  it  is  necessary  to  use  the  standard  engineering  driving 
force  equation  given  by  Sherwood  and  Pigford  (Ref.  16): 

w  -  (k4irr2)  (a  Pa)  (IV-49) 

Ranz  and  Marshall  (Ref.  11)  determined  the  mass  transfer 

coefficient,  k,  for  vaporizing  drops,  in  terms  of  the  Nusselt 

number  for  mass  transfer,  Nu  • 

w  ro 


( I  V- 5  0) 


(IV- 51) 


„  D  M, 

w  -  2nr  - L  Num  <*Pa  ( I V- 5 2 ) 

R  T 

The  film  thickness  may  be  empirically  determined  from 
equations  ( I V- 4 8  and  52)  as 

I1  -  Num  -2  ( I V - 5 3 ) 

Substitution  of  ( I V - S 1 )  into  (IV-52)  for  a  concentration  of 
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Substituting  equation  (IV-50)  into  (IV-49)  gives, 
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FIGURE  19*  Schematic  diagram  of  heat  tranafer  to  vapor 
film  and  liquid  drop.  (Reference  1) 
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This  empirical  expression  shows  the  dependence  of  the 
vaporization  rate  upon  physical  properties  and  upon  flow 
p arame  te rs  . 


Substituting  for  the  Reynolds  number,  the  vaporization 
rate  equation  takes  the  form: 
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If  the  liquid  temperature  and  vapor  pressure  do  not 
vary  with  time,  the  following  nondimensional  vaporization 
rate  is  obtained: 
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Equation  (IV-S6)  may  be  written  in  terms  of  the  Reynolds 
number  of  the  drop  based  on  the  speed  of  sound  (Red) . 
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Before  the  nonlinear  partial  differential  equations  can 
be  numerically  integrated,  the  boundary  conditions  must  be 
defined.  Since  the  response  of  the  system  is  required,  i.e., 
whether  the  system  will  amplify  or  attenuate,  the  conditions 
define  an  initial  value  problem.  The  boundary  condition  or 
initial  value  is  analogous  to  the  disturbance  induced  in  a 
rocket  engine.  Presently,  two  types  of  disturbances  arc 
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The  second  type  of  disturbance  consists  of  an  instantan¬ 
eously  imposed  velocity  component  in  the  ft-directim. 


ve  -  Ay  S IN 6  ( I V-63) 

At  constant  pressure,  temperature,  and  density. 

The  first  type  of  disturbance  is  a  pressure  disturbance 
and  corresponds  to  an  ideal  pulse  gun  disturbance.  It  should 
be  noted  that  it  is  possible  to  introduce  any  wave  shape 
desired.  The  second  type  of  disturbance  is  a  velocity  change 
and  corresponds  to  a  standing  or  spinning  transverse  wave. 

Work  done  at  Princeton  (Ref.  17)  shows  that  this  type  of  field 
is  produced  by  a  real  pulse  gun  disturbance  as  well  as  a  pres¬ 
sure  gradient.  A  third  type  disturbance  which  can  be  intro¬ 
duced  in  a  vortex  or  a  constant  tangential  velocity  disturbance. 
This  type  of  wave  most  closely  represents  a  gas  flow  disturbance. 
This  nonlinear  model  which  relates  instability  to  the  disturbance 
character  is  very  flexible  and  can  be  used  to  determine  the 
effects  of  shape,  amplitude  and  various  combinations  of  pressure 
and  velocity  disturbances. 

For  times  greater  than  t'-O,  the  pressure,  the  temperature, 
and  the  velocity  at  each  position  in  the  combustor  or  annulus 
are  computed  from  Equations  (IV-40)  to  (IV-47)  and  (IV-57), 

From  these  results,  stability  limits  can  be  computed  as  shown 
in  Figure  6. 

The  nonlinear  differential  equations  are  solved  numeri¬ 
cally  using  a  marching  technique.  The  computer  program  is 
outlined  in  Appendix  II.  This  numerical  technique  involves 
solving  for  derivatives  in  the  direction  of  march  (time 
direction)  in  terms  of  the  dependent  variables  (P',  T* ,  v' ,  p  * ) 
and  their  derivatives.  Derivatives,  with  respect  to  other 
independent  variables  (6),  are  obtained  from  values  of  the 
dependent  variables  by  a  finite-difference  scheme.  Using  the 
computed  derivatives  in  the  t  direction,  values  of  the  depen¬ 
dent  variables  may  be  computed  for  the  incremented  time.  The 
process  is  repeated,  and  dependent  variables  are  computed  for 
successive  values  of  time.  Initial  conditions  are  introduced 
in  the  form  of  a  disturbance  at  tine  ■  0,  By  analysing  the 
magnitude  and  rate  of  growth  or  decay,  stability  can  be  deter¬ 
mined. 

Stability  limit  curves  can  be  computed  as  a  function  of 
the  significant  instability  parameters  Av,  L,  and  Re^.  Re¬ 
sults  of  the  calculation  of  the  disturbance  amplitude  required 
to  trigger  instability  show  that  the  annulus  is  most  stable 
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when  Av  approaches  infinity,  Re^  approaches  zero  anti 
L  is  less  than  0.1  and  greater  than  5.0.  Second  order 
effects  have  been  shown  for  the  viscous-dissipation  parameter, 
specific  heat  ratio  and  axial  gas  Mach  number. 

Work  is  presently  in  progress  at  Dynamic  Science 
Corporation  to  include  the  effects  of  droplet  damping  and 
drift  terms.  Since  the  simplified  model  presented  here 
requires  the  use  of  an  instability  drop  average,  cevelopment 
of  a  model  based  on  po ly di s pe rs e d  sprays  and  bipropcllant 
droplet  sprays  is  under  development.  As  the  models  become 
more  complicated  simplified  nondimensiona 1  groups  become 
harder  to  use,  since  the  number  of  family  curves  required 
become  prohibitive.  Work  is  in  progress  to  integrate  the 
steady-state  model  and  the  instability  model  directly  so 
that  stability  zones  for  a  given  engine  may  be  generated 
by  computer. 


i 
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APPENDIX  I 

Dynamic  Science  Corporation  Steady-St att 
Spray  Combustion  Model  with  Two  Flame  Fuel 
Burning  for  MMH-NTO  System  Program. 


APPLnJDIX  I 


STEADY -STATE  PROGRAM 


1 .  Introduction. 

A  FORTRAN  IV  computer  program  has  been  written  which 
integrates  the  solution  of  the  simultaneous  first  order 
differential  equations  describing  steady -state  droplet 
vaporization  with  vapor  phase  decomposition.  The  equations 
involve  simultaneous  heat,  mass,  and  momentum  transfer  coupled 
with  varying  gas  dynamics.  The  equations  arc  solved  simul¬ 
taneously  at  distance  (x)  increments  along  the  axis  of  the 
chamber  for  each  of  a  finite  number  of  drop  groups  describing 
the  spray  distribution  of  the  propellants.  Paragraph  2 
describes  the  required  input  data.  Paragraph  3,  describes 
the  program  logic,  important  variable.',  and  numerical  inte¬ 
gration  technique.  A  program  listing  is  given  in  Paragraph 
4  and  the  solution  of  a  sample  problem  is  shown  in  Paragraph  5. 

2 .  Input  Data. 

The  computer  program  as  formulated  is  general  in  that  any 
propellant  combination  whose  burning  can  be  described  by  the 
model  can  be  analyzed  providing  that  the  thermodynamic  pro¬ 
perties  are  available.  The  input  data  to  the  program  includes: 

(a)  Propellant  thermodynamic  and  physical  properties. 

(b)  Chamber  geometry. 

(c)  Injection  parameters. 

Each  of  the  above  areas  will  be  discussed  separately.  The 
propellant  thermodynamic  and  physical  propellants  required  are 
tabulated  in  Table  I.  These  properties  are  input  as  function 
subprograms  and  are  curve-fit  to  yield  continuous  values. 

Table  I  also  shows  the  subprogram  of  the  respective  properties. 
The  program  listing  shows  the  property  subprograms  as  used  in 
the  analysis  of  the  MM11-NTO  propellant  system. 

Chamber  geometry  is  input  in  two  function  subprograms 
A(x)  and  AP(x).  Function  A(x)  describes  the  variation  of 
chamber  area  with  distance  (x)  while  AP(x)  specifies  the 
derivative  of  the  chamber  area  with  respect  to  distance  (x). 
These  functions  as  shown  in  the  listing  are  for  a  conical 
chamber  configuration. 

The  injection  parameters  are  input  in  card  form.  The 
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TABLE  I 


PROPELLANT  PROPERTIES 


NAME 

SYMBOL 

UNITS 

SUBROUTINE 

Liquid  Density 

Pi 

#/in3 

RHO 

Liquid  Specific 
Heat 

Cpi 

BTU/ #-  *F 

CP  L 

Vapor  Pressure 

PA 

#/in2 

PVAP 

Vapor  Viscosity 

wa 

#/in-sec 

VISCV 

Vapor  Specific 

Heat 

Cpa 

BTU/#- *F 

CVAP 

Vapor  Thermal 
Conductivity 

ka 

BTU/ sec-  *F- 
In 

FKA 

Heat  of  Vaporiza- 

X 

BTU/# 

tion 

PRO 

JUCT  PROPERTIES 

Adiabatic  Flame 
Temperature 

Tg 

*R 

TGINT 

Viscosity,  Thermal 
Cond.  $  Specific 

“*•  k*'  cp, 

VKSGAS 

Heat 

Dif  fus i vi ty 

D 

In2/sec 

DIFFU  ; 

input  variables  and  format  with  their  respective  units  are 
represented  in  Figures  20  and  21.  Specifying  these  parameters 
as  card  input  easily  permits  parametric  studies  to  be  made 
for  a  particular  propellant  system  and  chamber  configuration. 

Thus  the  effects  of  changing  such  conditions  as  initial  drop 
temperature,  injection  velocity,  0/F  ratio,  and  droplet 
distribution  may  conveniently  be  determined. 

3 .  Program  Nomenclature. 

The  logical  flow  diagram  of  the  steady-state  program  is 
shown  in  Figure  22.  The  FORTRAN  variable  names  describing 
the  important  parameters  of  the  model  are  shown  in 
(Program  Nomenclature).  The  equations  solved  at  each  calculation 
box  of  the  flow  diagram  are  described  in  Section  III  of  the  text. 

The  program  uses  a  second  order  Runge-Kutta  integration 
technique  which  is  described  below. 

3 £  ■  f(x.  y) 

and.  iL  .  £  “-f(Wn). 

Then  the  slope  of  y(x)  midway  across  the  Ax  incremert  (a)  is 
approximated  by 

[£]/  f<*  *  *f.  '  ♦ 

this  gives 


and  calculating  y  (x) 

yn+l  m  yn  *  c Ay) a 
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A  standard  card  form,  .IBM  electro  888157,  is  available  for  punching  sowce  statements  from  this  form. 


FORTRAN  CODING  FORM 


SECOND  VAPOR  I  EAT  1  ON  LOOr 
1  DANTI  CAL  TO  MIST  II- 
C1PT  USES  DIRIVAT1VSS 
CALCULATED  AT  AI/j  FOR 
SECOND  HALT  OF  RUNGS- 
I  IUTTA  INTEGRATION 


MINT  VALUES 
CALCULATE 
STASILITT  PARA, 


STOM  VALUES 
FOR  GRAPHING 


j  FIGURE  22 


Program  Nomenclature* 


FORTRAN  SYMBOL 

UTA  A 

AX  from  Fen  A  (X)  A 

v  '  c 


AXMIN 

B 


B 

C* 


CD ( RE )  from  Fen  CD(RE)  Cd 

CP 

Fen  CPL 
Fen  CVAP 

CPA  from  Fen  CVAP 
CP B  from  Fen  VKSGAS 
CPG  from  Fen  VKSGAS 
CP  AV 

DFU  from  Sub  DIFFU  D 

Do 

E 

F 

FVAP2 

G  g 

I 

XL  from  Fen  XK  K 

k 


DESCRIPTION 


nozzle  contraction  ratio 

chamber  cross-sectional  area, 
function  of  x 

area  of  throat 

boundary  layer  film  thickness 

characteristic  exhaust 
velocity 

drag  coefficient 

heat  capacity 
--  of  liquid 
--  of  vapor 

--  of  vapor  at  average  mantle 
temperature 

--  of  bulk  products  at  average 
mantle  temperature 
--  of  bulk  products  at  static 
temperature 
--  of  mantle  mixture 

molecular  diffusion  coefficient 

injector  orifice  diameter 

activation  energy 

injected  fuel  fraction 

local  fuel  fraction  vaporized 

acceleration  of  gravity 

momentum  flux 

evaporation  rate  constant 

thermal  conductivity 


*Fcn  means  function  subprogram 
Sub  means  subroutine 
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FORTRAN 

FCN  FKA 

SYMBOL 

FKB  from  Sub  VKSGAS 

FKG  from  Sub  VKSGAS 

FKAV 

1 

EM 

a 

RVAP 

m 

WTML  from  Sub  WTMINT 
WTMAV 

WTMOL  (input) 

M 

DN 

n 

FNUH 

Nuh 

FNUM 

Nua 

0 

FVAP1 

e 

P 

p 

PA  from  Fen  PVAP 

p. 

PR 

Pr 

SDSPD 

a 

<lv 

RS 

r 

rmo 

DESCRIPTION 
--  of  vapor 

--  of  vapor  at  avara|«  aantle 
tarparaturt 

- -  of  bulk  products  at  avaraga 
■ant la  taaparature 
••  of  bulk  products  at  static  temp. 
- -  of  aantle  mixture  at  av.temp. 
length 

aass  of  drop 

vaporization  rata  of  spray 
(fraction/in.) 

■olecular  weight  of  vapor 
--  of  aantle  Mixture 
--  of  propellant 

nuabar  of  drops/sac.  in  each 
drop  size  group  or  reduc* 
duction  velocity  constant 

Nusselt  number,  heat  transfer 

Nusselt  number,  aass  transfer 

injected  oxidizer  fraction 

local  oxidizer  fraction  vaporized 

static  bulk  gas  pressure 

vapor  pressure 

Prandtl  number 

speed  of  sound 

heat  of  reaction 

heat  arriving  at  drop  surface 

drop  radius 

1/m 

tn  r*/in  *]  summation  over 

all  size  groups 
when  m  ■  2;  area  mean 

m  ■  3;  volume  mean 
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FORTRAN 

SYMBOL 

DESCRIPTION 

RM 

rm 

radius  oi  mass  median  droplet 

rmax 

maximum  droplet  diameter 

RL  from  Fen  RT 

rt 

distance  from  drop  Surface 
to  decomposition  flame 

R 

R 

universal  gas  constant 

RE 

REP 

REN 

Re 

drop  Reynolds  number 
--  bulk  gas  properties 
--  mantle  properties 
--  based  on  speed  of  sound 

SC 

Sc 

Schmidt  number 

TG 

TL 

TBR 

TSTG 

T 

temperature 
•-static  gas 
--  liquid 

--mean  mantle  temperature 
--stagnation 

TB  from  Fen  TBL 

Tb 

boiling  point  temperature 

TWB 

Twb 

wet  bulb  temperature 

U 

u 

gas  velocity 

VD 

V 

u  -  V 

ABSF  of  VD 

A  V 

|u  -  v| 

V 

V 

drop  velocity 

V(input) 

liquid  injection  velocity 

W 

w 

vaporization  rate  of  single 
drop  (lbm/sec) 

EMDOT  (input) 

• 

W 

flow  rate  (lbm/sec) 

X 

defined  by  equations  27  $  28 

X 

X 

axial  position  (x-0  at 
injector) 

z 

defined  by  equation  (correction 
factor  for  mass  transfer) 
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FORTRAN 


DESCRIPTION 


RH  from  Fen  RHO 
RHOG 

RiiOAV 

GAN 


Fen  ALM 

VIS  from  Sub  VKSGAS 
VISB 

VISAV 


RAT 

SUBSCRIPTS 


SYMBOL 


a  defined  by  equetion  (6) 

0  reeetlon  velocity  constant 

p  density 

- -  of  liquid 

- -  of  bulk  gas  at  static  tenp. 

--  of  aixture  in  mantle  at 
average  teaperature 
Y  ratio  of  specific  heats  of 

product  gas 

efficiency  based  on  characteristic 
velocity 

X  heat  of  vaporization 

v  viscosity 

--  of  bulk  products  at  static 
teaperature 

- -  of  bulk  products  at  average 
mantle  temperature 
--  of  aantle  nix. at  aver. temp, 
v  kineaatic  viscosity 

0  surface  tension 

e  aixture  ratio,  O/F 


a 

c 

d 

f 

2 

i 

i 


0 

p 

v 

6« 


vapor 
combusted 
drop  let 
fuel 

product  gas 
at  injector 
liquid 

oxidizer  or  stagnation 

particle 

vaporised 


n  n  n 


Program  Ustinj 


PROGK AM  VMn 

STEADY  s T « I  l  Lui>iuubl  iuiv  l*iuutL  —  T  rtu  FLaMl  i-ull  dukNIimo 

DYNAMIC  SClENCL  LOKKUrtHfiu^ 
l90o  toALfvL*  AVl.  NUb 
MuNkuVIa,  i  mL  i  r  JKIX  1  M 
0  5  7  -  Ji5i  (M.t.  2  1  J  ) 

OK  tt.N  sObVlx. 


MuNkuVJa,  i  mL  i  r  JKIX  i  M 
0  5  7  -  3  C3  X  (  M  •  v.  •  2  1  J  ) 

oKLlN  sObVli. 

COMMON  OR 

O  1  M  t  Nb  i  u  x  Hulijy),  uKli .  5  u  x  >  u 

l-  1  Me. NS  i  ON  i.od  (  9  ) 


DC  i>X  l/)i  oCuY  (  b  ) 


TYPb  INTF.olk  hi.  dX  »  nv.  >jY  •  lm.<1  x  dijI  xdl>3 

0  I  Mb  NS  I  uN  T1TLl(3v)  ♦TL(4(J)»  V(40)x  lMOOT  (  2  )  ♦  KM  ( 2 )  x  b  *  u  (  c  )  »  Y  (  5  )  ♦ 

1  W  T  muL  (  2  )  ♦  mk  (  2  )  *  SC  (  2  )  ♦  NSb  T  (  2  )  »  F  A  (  2  )  *  LAN  I  (  2  )  *  t-  n  M  (  c  )  •  r  Mi  (  3 

2  tv  C  U  N  (  2  )  x  K  b  b  (  2  u  )  x  L  M  S(20)*bM(4j)x  P  l  (  <:  )  x 

3  I  4u  )  •  JN  ( <4 u  )  x  ipix b  (  2 0  )  x  t  NkiP  (  2  )  x  VO  (  4(J  )  x  Li*'  i  ix  1  (  <* 0  )  * 

4  VLlUOf(4U)»L'Vl40)»L'lvl(4o)»L>T(4ur»DLTA(4u)*br'.l£)»jrO(£)» 

P  bVft  U  )  tbxK  U  I  »  jui/K(^  )  >iD\/K  J2  (  ^  )  *okl\(  :  )  » 

<3  Kio(2)*Kbu(2)*rlVb4l(t)*KDVb^i(<;)*  Rivil  aix  <  O  )  x  k  L  ix  (  J  b  ) 

QUO  i  VAllNOL  (  KMt  aN  x  R  1 0  )  x  (  kMlaix  (  3  )  x  R3Q  )  x  C  K_Mt_AN  (  5  )  »kVjc  1  )  • 
i  “  ( RMlaN ( 7  j x  RdV  32 1 ) 

L)  I  Mb.  NS  I  ON  KM  20  >  xXbNTR  <  20  >  xWbNTRI  20  )  x  TwbINT  (  20  >  xObLXI  23  !  x 
1  OUMI TL( 2U xDiMTL ( 2CV 
01  Mb  NS  I  ON  VD10 ( 2 )  xVDpO  <  2  >  xV0321 ( 2 )  .VD13  32  1  (  2  ) 

DTWENblON  wM20) 

OlMtNSIuN  iu 1 Sx ( 40  I 


OAT  A  (  R- 18  540.  )  ,  (  TWOP  I=to.2631b33)x(G=3bb.4)x(FTPi=4.1bb7902) 

IN  Lbb/MOLt  R  IN/SbC/SbC 

"V'KtAD  (  5'»  10  )  (  T  I  TLt  (  I  »  Vl  =  1  xo9)' 

lu  hQRMAT ( 13Ab )  _ 

I F ( tOF x  3  )  i  3  x  1 6 
13__bToP _ _ 

16  WRI Tb (6.2u) ( TI  ILL (  i  )  t  f=i x39) 

2u _ FORM AT ( 45H1 SN 70u2  S  T  t  AOY  — S  I.aT  t  SPRAY  ljMdUSTI JN  i»iuolL  /  (  j.  h  IpAb 
RbA0(  S»30)P.RAT  .CAM,  HL(I)*  V  (  I  )  ♦  bMOOT  (  I  )  ,Kh  1  1)  .  b  I  O  (  i  )  ♦ 

1 _ w  T MQL  (  1  >  x _ NSbT  (  I  )  x  I  =  1 »  2J _ 

30  FORMAT (3F9»L/!6F9xCi 14) ) 

WRJ.TfcA6»4.U)  (bMOOT  (  1  )  *  RM  (  I  )  ♦  S I  o  (  I  )  *  WTMOL  (  I  )  x  1  -  1  ♦  ) 

40  FORMAT ( 1H0 

— i— §21-  &0H  MASS  FLOW  __R_M_.. _ Sj.G  MOLL  WT  / 

1 2  FI  O  FI p . 2  x  L 1 2  •  3  x  cF7.2  / 

_ £ _ LULa-Li _ LUL1A1 _ ZL1 ,2 _ )  

READ  (S»41)  RATMN.RATMXxAXMlNxAXMAXxOFSTOC 

41  .....  f.ORM A.T  J  3  F  1  p .  3  )  _  ..  .  . 


0  =  0  • 

lbZ  =  u - 

J  L  =  1 
X"=0'i  ' 

FA ( 1  )=0. 
"TATZT^TTi" 
PHI =0. 


.  =  OAM-J.x 


(jAM2  =  GAM1/2. 

"  ~G  A"m  Y=  G  A’ M  2"7  oa  m~ 
bAM4=GAMl/0AM 


GAM5  *  GAM/GAMI 
'  r«D'*EWD'Om  f+EMUOTm 
THRD*l#/3. 

- ITT-I -  , 

KSW*  1 

C  . . - - - - - 

c 

DO'TUO'TilVZ"-" - - - -  ‘ 

DNT ( I ) =0* 

WdON ( I ) =TWOP I *WTMOL ( I ) /R 
N*NSET  (  I  ) 

CTL"l_"l5SnT?hTRm  TT)' VsTg'iT)-#  n7  r'ssT 

CALL  MASSET(EMS*RHO(TL<  I  ) . I ) . RSS . tMDOT ( 1 )  •  N ♦ DNS ) 

- N5FnTT'i‘_N+4 - 

N  =  NSFT( I ) 

- DO  100  J=T •  hi - 

K*  J  + J-2  +  I 

V<  K  ) =V ( I ) 

- en"T<r"T=rMSTjr - 

EMINT(K)  *  EM(K) 

- R?r  g"T=RSSU) - 

ON (  K.  ) =DNS ( J ) 

o"NTrTr*_oNTrrr+ONfsTjT 

KK ( K  )  =  1 

— . FwTKT"=""oV - 


100 

CONTINUE 

c 

ONTT*DNT< 1 )+DNT(2 ) 

XT  =  5  • 

DX  =  6. 

RV'AP  *  G.20 

GO  TO  106 

C 

C 

*****  END  INITIALIZATION  -  START  CALCULATION 

C 

C 

#*********#**************#****#**##*###*#*###*******#*#**#*##.###*# 

C 

C 

********#****«*******#**«***«##**##*****#**##*#####****#«  *#**#«»*# 

"TW 

105 

lF(X+DX-XT)105. 105*999 

IFIX-O. 25)110*110*107 

106 

50 

READ(5.50)DX.XT.IP,IPLOT 

FORMAT (  2F9*  0  *219) 

107 

710  fO  108 

OX  =  XDEL(RVAP) 

ToF 

0X2  *  DX/2. 

IF  ( DX ) 999*999* 110 

110 

P2=P*2 • 

P2R=P2/R 

c 

c 

*****  OBTAIN  STAGNATION  TEMP.TSTG.  ANO  THt  DtR I  VAT  I Vt  OF  TtMP 

c 

WITH  respEct  to  o/f*  tgod 

CALL  TGI  NT  ( RAT *P . TSTG . TGOU ) 

-c- 

c 

*****  OBTAIN  MOLECULAR  wti&HT  OF  PRODUCTS  AS  FUNCTION  O/F 

c 

CALL  WTMINT  < RAT *P . WTML . WTMLDT > 

rg=r/wtml 

U2=U*U  *• 


’C 


*****  CALCULATE  STATIC  TtMP.  FROM  TOTAL  ThMP.  AND  GAMMA 


...j 


ro=TSTO*(  l.-0AM3*U2/KG/TSTG/u) 
C 
C 

c  *#**#  CALCULATE  GAS  DENSITY 

RHOG=P/KG/TG 
C 


L 

c 

v. 

u 

c 

L 

C 

C 

C 


C 

C 


L 


C 


C 


*#*'**  CALCULATE  VISCOSITY*  TntRiMAL  CONDUCTIVITY*  ANu  sPcCIFlC  nt 
OF  oAS  AS  FUNCTION  uF  P,  TG*  O/F 

'  “'Call  '  vksJas  '(  ra7 »p,TTT*vi  sVFku*cpg ) 

kho'VS'=RnOo/V  1  S*2  . 

KHOT  E  =  •  3  7  P  ^KHUu 


***'»****»il*:»lHrli 


*U****iHHHHhtinHM**t***** 


*****  STAkT  DROP  VAPORIZATION  LOUP 

12o  DO  200  1=1,2 
N  =  iMSET(  I  ) 

ENMP(I)=C. 

10(3  2  00  J=  i  » N 
K=J+J-2+I 

<|3TSC"(K?  s  0 


*****  TEST  FOK  ZERO  MASS  OF  DROP  -  SKIP  IF  ZERO 

_ I  F  <  EMI  K  )  )  _  2  0  0  ,_2u  0  ,  i  3  U  _ 

T 3  u  _  T B  R  =  (  f  G+T~l1~k17/2  .  . 

_ C  A L L  V_K S_G  AS  ( l<  A  T  iP  »  T oK  ,  V_I  So  ,  F_K. b ,  C P u  ) 

*****  qotIAN  \ApQKlZATluN  PKtbSUKE  ,PA,  AT  LluUlu  T  •  *  1  L  » 

PA  =  PVAP  (  TL  (  K  )  i  1  )  '  ~ . 

_ _P  A  2_P=_P  A/ P_2 _ 

” ’dp  A~2P"=l"'lp  A 2  P 

_ C_PA_=C_VAPJ_IBR  ,_I  )  .  ___  . . . 

RH=’rh6(TLK7,T) 

_ R V  =  W TN1UL  (  I  )  »P/R/TL  (  <  ) _ 

VD(K)=U-V(<) 

RS(K)=(EMIK)/FTPI /Rm)**ThkD 


_C _ t YNOLD_S..NUMBcR  OF  DROP  AT  UAS  PROPERTIES 

""RE*  RS~(  k1~*Rh67s*ASS~F  (' v'l5'(V)  ) 

C _ 

to TMAV  =  PA2P*WTMOL( 1 )  +  DPA2P*WTML 

_ VJSAV  =_  VIbCy|_T_bR_,  n_*_PA2P_  +  DP A2P*  Vj.S.b _ 

Fk"av"-='"PA2"p#_FKA  (YbR~»  I  r+  DPA2P~*FKB 

__RHOAV  =  P*WTMA_V/K_/J_oK _ 

CPAv'=  '{PA2P*WTMOL(  I  )  *CPA~T~DPaYp*W  frtL*CPB  )  /  to  TMAV 
REP=  2,*RS(  K  I  *RFiOAV*AbSF  (VD(K)  I/VISAV 
PAIR  =  PA  /  2. 

CALL  DIFFU  <  R  A T , P , TbR , PA  I P , I , DFU ) 

- DwvDTuTfpyrbP - 

PR(I)  =  CP AV*V  I  SAV/FK AV 
5rTTr=-vrsAV7ffnuAv7u-Fu 
RTRE  =  SQRTFIREP) 


FNUM  =  2. 

+  0.6  *  SC  (  I  )  **THRl) 

»  RTRE 

FNUH  =  2. 

+  0.6  *  PR  (  I  )  **TFiRD 

*  RTRE 

- _c - 

c 

69 

4  y?r~**e 


Z - CHECK  FOR  ENTRANCE  TO  blSSOcTZnTOTi  LO UP 

GO  TO  ( 149*160)  I 
16U  B=  2 «  * R S  <  K ) / ( FNuR-2 • ) 

IF  DIVIDE  CHECK  161*162 

T6T'H‘V-i"F3fc”‘” -  - - 

162  KL  =  RT( RS(K)  ) _ 

IT  (  B-R L  )  149*  lb3  »  16j 


L  ENTER  DISSOCIATION  FLAME  LOOP 

~  T61  ~  T  F  kk  "f k  -  - 

KOISCIK)  =  1 


00  TO"  (  164,1661  .  I  f 

-  - 

164 

KK ( K)=2 

c 

iTJHilVAnnniKiftifnnnifltMPiViTniBimn/nHnnin 

I F ( X  )  601  *601,600 

| 

6  OTJ 

XENTR(K)  =  X 

WtNTR(K)  =  wW(K) 

XL  =  XK  f  R S ( K  )  ) 

W  =  XL  *  TWOPI  *  RS ( K )  *  RHO ( T L ( K ) , I  ) 

Ll  =  W  *  CPA  *  RL  /  FKAV  /  2.  /  TWOPI 

GO  TO  602 

/  RS  (  K  ) 

/ 

RSI  <) 

601  XENTR(K)  =  -.001 

XL  =  XK ( R S ( K  )  ) 

WtNTRT  K )  =XT>  TWGFI  *  KSCikF  *  KnoTTLfK)»TF 

LL  =  WENTRU)  *  CHm  *  Kl  '  KAV  /  2.  /  TaoPI 

/  Kb ( K  ) 

/ 

K  S  (  K  ) 

“5TJT 

TWb I  NT ( K )  =  BTTTu 

DELX(K)  =  V(K)*EM<K)*CPl(Tl(K)  *1  )*<  TWblNTU) 

—  T  L  (  K  )  )  / 

1  TWOPI  /RS ( K )  /RHO ( TL ( K ) , I )  /  562.  /XL 

0UM1TLIK)  =  TL ( K ) 

OUM  T  L  (  K  )  =  T  L  (  K  ) 

MAIN  DISSUOiATION  1- LAME  CALCULATION  LOOP 

16b 

XL  =  XK(RSIK) ) 

W  =  XL  *  TwOPI  *  Kb(K)  *  KHu ( T L ( K ) » I  ) 

1 

LL  ~  W  *  CPA  *  KL  /  FKAV  /  2.  /  TWOPI 

TWB  =  817.0 

/  KS(K) 

/ 

Kb  I  K  ) 

c 

TEST  FOR  X  GREATER  THAN  XENTR  +  LENGTH  TO  REACH  WET 

bULb 

XTRA  =  XENTR(K)  +  DlLX(k) 

1 

IF  (  X-XTRA ) 166, 167 , 16/ 

166 

CONTINUE 

1 

-  -.1 

c 

I N T C.RPULA T c.  To  k1imi>  AO  Ainu  Two  wHlN  X  L.ln. 

XliUK  + 

UbLX 

W  =  WENTRIM  +  <W  -  a  l.  N  T  K  (  K  )  )  /  DELX(K)  *  (  X-XENTk  I  k. )  ) 

! 

.....  -A 

TWb  =  DUM1Tl”(K)  +  (  Tao'-OOKii  TL(k)  ) /DeUx  (K>*(  X-XENTK  (K  j  ) 

OT ( K  )  = i T WO-OUMTL ( K )  1/0X2 

DUMTL(K)  =  TWB 

DM ( K  )  =  -W  /V<K) 

c 

> 

BACK  TO  MAIN 

BETA  fk)  =  W  /  R 6  (  K  )  /  FNUM 

GO  TO  156 

* 

c 

CALCULATIONS  WHEN  X  GR.TH.  XENTR  +  UELX 

16/ 

DM ( K  )  =  -W  / V ( K )  T- 

DT  (  K  )  "  (  T  WB-DUMT  L  (  K  )  )  /  U'X  2 

DUMTL(K)  =  TWB 

bETA(K)  =  W  /  RS ( K )  /  FNUM 

r 

GC  TO  156 

149  PMDP=0.99*P 


13^' 

IF (PA.ot.PMOP )  130*  153 

Td  =  TbL ( P* I ) 

WKI  T  t ( 6  »  2  OOU  ) 

7  GOO 

C 

FOKPlM  1  (  £  lrtuG  TA  1  lMci*  T  wUi'idcK  130) 

U  =  0  ,| 

W  1  =-i)M  (  K)  *V  (  V  ) 

Wl  =  ABSF (Ml) 

RRS=Rv*R3  (  k  )  *RS  (  k  ) 

cUE"=Tnum»f  k. 4V#RS  (  K  )  #T  wOP  1  /CPA 

COc!  =  AbSF(COL) 

152 

W= C0E*L06F { 1 .+CPA* ( TG-Tb ) / ( ALM<  Tb» I )+8. 76E-10* (Wl/KRb ) **2 )  ) 

C 

WR  =  ABSF  (  ( Vw  —  W 1 )  /  Wv  ) 

15  3 

XT  =“  K 1  +T  "" 

I F ( K1.GT.15) 154,133 

w  I  =  f  w  +  w  f )  72 .  o 

IF  ( wR.Lt.0.02 )  134, 1 3  ^ 

1  3  4 

DM  (  K.  )  =-W/V  (  K  ) 

D  T  (  K  )  =  U  »  0 

TL(iO=Tb 

Dt  T  A  (  K. )  =  W  /  R b  (  K.  )  /  FNUM 

135 

GO  TO  136 

CCNT I NJt 

o  - \ 

1 

• 

CALL  COMP  (KAT.P.Y) 

PL  (  1  )  =  0.0 

34  OU 

GO  T  0( 3400,3410 IKbW 

PL ( 1 )  =  O.C 

341  U 

PE ( 2 )  =  0.0 

b E T A ( K )  =  WCON(  I  )*DPOT*LOGF( (P-PE( I  )  ) /(P-PA) ) 

C 

W=6ETA(  K.  )  *RS  (  K  ) *F  NUM 

C 

1 F  ( W )  190,  1  9  ^  »  192 

19U 

W  =  0. 

WRI TE(6,2050 ) 

2u5u 

C 

FORMAT(19HOW  NEG.  SET  TO  ZERO) 

C 

DM  (0=0. 

QT  (  K  )  =Tl*uP  I*FNUrl*M<l«V*(  To-TL(  K)  )  /  Ert  (  K.  )  /CPL  (  TL  (  K  )  » 1  )  /  V(  M  *Kb(  lU 

GO  TO  156 

C 

r 

V. 

c 

192 

DM ( K  )=-W/V(K ) 

Z=W*CPA/TWOP I / RS ( K  ) /FNUh/FKAV 

EZ  =  EXPF  ( Z  )  ~1  • 

• 

1)1  F  =ALM  (  TL  (  K  )  ♦  1  1-CPA*  (  l^-TL(K)  )  /tZ+8#78t-10*  (  w/RV/KST  K  )  /RS  (  ^T)**2 

DT  (  K  )  =-D  I  F  /CPL(TL(<)  » I ) /EM ( K ) *W/V ( K ) 

c 

• 

IF( DT(K) +100. ) 7711 , 156,136 

77TT  UTTRT  =  -  IXTC  . 

c  i 

c 

c 

1 

1 

7; 

-  - r5-g. 

— caRTWor  .  . 

DV(K  isRH0TE«VD(O*AdSF(V0(K.)  )  *CD  <  RE  ) /RFi/V  (  K  )  /  RS  (  K. ) 

n\  r 


ENMP ( I  )=ENMP( I )+DNU)*DM(K) 

WWVRY-*  W 

CONTINUE 


KK  =  RG  *  T 0 

505PD  =  SuRTF  <  .^K  *  oAM  *  o) 
1 F (U.GT.SOSPu)2ub»21(; 


NP  =  0 
AX  =  A (XI 

PrtIP  =  -  (  lNmP  (  1  )  -n.i*i*ip  ( 2 )  ) 
PflT=FATIT-FFA7  2  ) 

U  =  PHl/RhOG/MX 


EMACH2  =  tMACn  *  #2 
UO  TT7  7  23u »  2  20  ) 

RA  TP  =  (  F  A  (  1  ) /F  A  ( )  *tNMP  (  2  )-f  MP(  1 )  ) /FA(  2  ) 
IFfRAT. GT  VKATMN.ANC  •RAT.  L  T  7k  AT  M  X  )  72~4T "727 

RATP=0 • 


GO  TO  240 
7CSW=2 
NP  =  IP 


Hr  #  #  *  *  #  ■! 


Tf  7"NP~_  -  “  i'2" 5>  0 . 3  0  0  *  3  0  0 


If#*-**-##"****  *  #***'#»•«■**• 


'Continue 

NP  =  NP+ 1 
X  =  X  +  DX  2 


00  2b 'J  1  =  i.  *  2 
T  A  (  1  )  =  t  M  0  o  T  (  i  ) 
N  =  Nb  t  T (  1  ) 


00  260  J= 1 ♦ N 
K= J+J-2+ I 

IF  (  t  M  (  K  )  )  2b  C  *2b0  *<;35 


2bb  tM(K  )  =  t  M  (  K  )  +  0  X  2  *  0  v.  (  K, ) 

F  A  (  F)  =F  A  (  I  )  -ON  (  K  )  * t M  l  <.  ) 

TL  (  K  )  =  T L  {  K  )  +  0X  2  *0  T  (  x.  ) 

V ( K ) = V ( < )+DX2*DV( < ) 

C  - 

C 

I  F  (  E  M  (  K  )  )lbu»ltiU»2bO 
lbo  cM ( K. ) =6. 

TL ( K ) =0 • 

v(o=o. 

RS  (  K. ) =0. 

Vl)(  K  )  =0. 

ON T (  I  ) =0N T ( I  ) -DN ( < )  72 

DN  (  X~)  =0.  " 

26U  CONTINUE 


KAT=FAri  )7FA  I  2  ) 

IFIRAT.LT.  RATMN)  261.  262 


GO  TO  264 

'COTCTIMDE  - 

IF(RAT.GT.RaTMX)  263.  264 

-T?AT‘*="RATMX 
CONTlNuc 


P2=P*2» 

P2R=P27R“"'  -  - - - 

CALL  TolNT  ( K A T ♦ P » T L i G . I GuD ) 

- XATL^wTFrr  i\T  ■  TR'AT  VPVwTPi  L“»  wT  ML  DTT 

i 


U2=U*U 

'TGVTFTG"*rrV-GWJ*'a27'KGVTSTG7“G') 

RHOG=P/RG/TG 


CALL  VKSuAS  (RAT .P.TG.VIb.FKG.CPG) 

"R  H0V3'=PW0g7  \7T  G'F"2  - 

KHUT  t=.376*KHOG 


SUMP  =0 

30^p-A"V--0-.-0- 
SUMPB  =  0.C 


WDSUM  =  C.O 


♦a************************************  **************************** 

112U  UU1200  1  =  1. ^ 

N=Nsrrrrr .  . 

tNMP I  I ) =0. 


DO  1 2 00  J=1,N 
__  K  =  J  +  J  -  2  +  I 

kdTscTk”  o' 

I F ( EM ( K ) ) 1200. 1200. 1  130 


1130  TBR= ( TG+TL ( K ) )/2. 

CALL  VKSGAS  (RAT  .P.TBR.VISb.FKd.CPB) 


PA  =  P VAP ( TL ( K )  .  I  ) 

_ PA2P=PA/P2 

DPA2P= 1 • -PA2 P 
CPA  =  CVAP( TbR, I  ) 


RH  =  RHO(TL(K) , I  ) 

R V=WTMOL ( 1 )*P/R/TL(M 


VD(K)=U-V(K) 

RSI K )  =  ( LM( K) /FTP I /RH  )**ThRD 


RE  =  RS(K)*RHOVS*AdSF(VD(M ) 

WTMAV  =  PA2P*WTMOL( I )  +  DPA2P*WTML 


- V1SAV  =  \71 51V(  TBR .  I ) *PA2P  +  DPA2P»VT'5b - 

FKAV  =  PA2P*FKA ( TBR. I )  +  DPA2P*FKB 

RKO'SV  =  TJ*'WTWAV?‘K77BTF 

CP AV  =  ( P  A2P*W  T  MOL ( I ) *CP A  +  DPA2P*wTML*CPti ) /WTMAV 

i 

1 

I 

- RIT>v--2-;ifR-5-ricr*'PFTUAV»7TB5Fll7i;rri;TTyVT5'AV' . 

PA I P  =  PA  /  2. 

MJia. 


DPOT  =  DFU*P/TBR  7, 
"P’PTTT"=""CP*AV»vT5WFTirAT' 
SC  C I )  =  VISAV/RHOAV/DFU 


RTRE  =  ^QRTf-  (REP  ) 

FNUM  =  2.  +  0.6  *  EC  <  1  )  «  *TmRl.  *  RTkl 
FNUH  =  2.  +  0.6  *  PKtl)**TMRD  *  RTRt 

C 

C  - 

C  CHECX  FOR  ENTRANct  TO  uISSOCIATION  LOOP 

uO  TO  (1  14* .1160)1 

1160  6=  2.  *RS(IO/(PNUfW.  ) 

IF' DIVIDE  CntCK  ll6l7lio2 

1161  d  =  l.tiO 

1162  KL  =  RT ( RS( X ) ) 

IF  (B-RL)  1  14*. 1 16j . 1 16 j 
C  ENTER  DISSOCIATION  HAMt  LOOP 

1163  II  =  KK ( K  > 

'XDTSCIX)  =  1 

UO  TO  (1164.1x03)11 
1  164  XX  (  X  )'=2 

o  CALCOLATt  mNu  STukt  IimITiAl  uUAN  1  I  T  1  c  6  oN  1ST  L  Y  L  L  L 

IF  OT)  611.611.610 
6 1  u  XENTR(X)  =  X 

WEN  !  R  <  k  )  =  Ww  (  x ) 

XL  =  X  K  (  R  S  (  X  )  ) 

W"=”XL  *  T  WOP  I  *  K6(X)  *  R liu  (  T  L  (  X  )  ♦  1  ) 

LL  =  wLNTK(X)  *  CPA  »  Ki_  /  FxmV  /  l*  /  T  a  jP  1  /  kMx 

oo  7 a  viz 

611  XtNTR(X)  =  -.001 
XL  =  XX ( RS ( X ) ) 

W,ENTR  (  X)  =XL*  TWOPl  *  koUI  *  KHu  (  TL  (  X  ).  I  ) 

LL  ='  WtNTR ( X )  *  CPA  «  kl  /  FXAV  /  4.  /  TWuPI  /  KS ( X 

612  TWblNT(K)  =  817.0 

"  DtLX  r<T  V'V  (  X)  *tM  (  X  )  *CPL  <  TL  <  X  )  .  I  )  *<  Ubl  N  T  (  X  )  -Tl  lx)) 
1  T  WOP  1  /RSI  X  )  /Ri-m<  TL  <  x  )  .  1  )  /  564.  /XL 
uUMlTL(X)  =  T  L ( X ) 

OO.’iT  L  (  X  )  =  T  L  (  X  ) 

x  MAIN  DISSOCIATION  c  LA, Ml-  CALCULATION  Loop 

1166  XL  =  X<(RS(X) ) 

W  -  XL  *  TaOPI  *  R.'(x)  «  Rho(TLlX).I) 

LL  -  W  *  CPA  *  ix  L  /  c  X  Ay  /  2 .  /  T  P I  _/  K  5  (  x 

Twa  =  H’17.0 

x  TcST  FOR  X  oREA  T Lk  I-iAk  a  l  N  I  k  -  llNuTm  Tu  klAlm  .vlT 

XTRA  =  XENTR(X)  +  ^lLX(x) 

IF  (  X-XTRA ) 1166. 1 167 , i 167 

1166  CONTINUE 

C  INTERPOLATE  TO  FIN  aO  m,ND  Two  amlN  X  L.Tii.  X  ENT  R  + 

W  =  WENTR(x)  +(W  -aENTR(K)  )7  DEl.’X(V)  *  (  X-XZNTR  (  X  )  ) 
Two  =  DUV I TL ( K ) +  (  TWc-OU Mi TL ( X  )  )  / utLX ( X ) *  ( X-XlNTk l X) 
D  T  (  X  )  =  (  Tlr,o-bUMTL  (X  )  )'/DX2  -Dt(.<)/2. 

DOMTL(K)  =  TWb 

DM(X)  =  -w  /Vi  X)  -  DMiO/2. 

DMC  =  -W  /  V(X) 

C  oTC^T  Tu-  M  A 1  N 

oLTrt(X)  =  V,  /  ka(X)  /  FNuM 
GO  TO  1156 

C  CALCULATIONS  WHEN  X  GR.Th.  XENTR  +  i/ELX 

1167  DM ( <  )  =  — W  / V ( X )  -DM ( X ) / 2 • 

DMC  =  -W  /  V ( x ) _ 

DTIKT  =  (Twa-DoMTL fx ) )  /  0*2  -07 ( X ) /l. 

DUMTL(X)  =  Twb  74 

BE  T  A  (  <  )  =  W  /  Rs  (  X  )  /  FNUM4 
GO  TO  1156 


/  IN  6  (  X  ) 


/  R  a  (  X  ) 


_^/_  K  o  (  X  ) 
aOLa 

ULLX 


c 


c 


1149  PM3P=0.99*P 

1  f  TP  AT  Gc  .  P  Ml)  P") "  T I  5  0",  'Tib'S 

1150  T B  =  TBL ( P  »  I  ) 

•  -  -icr'="o -  ""  - - - 

_ W1=-DM(K<»V(K) _ 

Wl  =  AbsF  (  w  l ) 

RRS  =  RV*KS( K ) *RS<  K ) 

'C0'£  =TNUri*  F1CAV*  R  ST  <T  #Tw  oPT  7  CPA" 

COE  =  ABSF(COE) 

'IT5T'"w='CU'EVCD"6rfr."+'CPTiTTcr-TeT7TAL«7TBVrr+‘tf;T8r-T0VrwT/HKb')**2n - 

WR  =  ABSF (  ( W-Wi  ) /to ) 

- ki  =  rrr  i - 

IFUl.GT. 15)1154, 1153 
"IT5T  wT  ='Tw  +  wTT/2.0 

I F ( WR.LE.0.02 ) 1154,1152 
"IT54"I5‘Rr"=_-wY"/"vTic7 

DM ( < )  =  DMC  -  DM ( < ) / 2  • 

DT ( <) =0.0 
TL ( K ) =TB 

'b'Et'aTk  )""=■"  W  "7  "R"S( '<)'"/""FN“UM'  . 

GO  TO  1156 

T - - - - - 

1155  CONTINUE 

- CALL  COMP  (KAT.PtYJ - 

PEI  1)  =  Y ( 5 )  *  P 

- - - -prfrr-r--C; - - - 

BETA(K)  =  WCONI I )*DPOT*LOGr I ( P-PE ( I ) ) / (P-PA ) ) 

- ¥-FETATrr_R---Y1>T7-_w - -  - - 

IF(W)  1190,  1192,  1192 _ 

1190  W  =  0. 

DM ( K ) =0. 

DM'C  VoTo"  .  . 

DT(K)=TWUPI*FNUrl*FKAV*(TG-TL(IC)  )/EM:m/CPL(TL(<)  ,  I  )/V(M*KMM 

"GO"TO  IT  5  6  . 

c 

1192  DMC — ^WTVTTj - 

DM(K)=DMC-DM(K)/2. 

-z-=-w'ic"PA/fWOPr/RST  K  ) /FnOh/FKAV 
EZ=EXPF(2 )-l. 


U I F  =  ALM ( T  L I K ) » I  )  -CP A* ( TG-TL I <)  ) / t Z+8 . 78 E- 1 0* I W/K V / Kb ( * ) / KS l n )  )**2 
DT  (  K  )  =-D  I  F  /CPL(TLIK)  *1  )/EMHC)*W/V(M  -DTIO/2.  _ 


- Z - 

IF(DT(K)+100. )7712* 1156, 1156  . 

7712 

C 

DTIK)  =  -  100. 

C 

1156 

**«****##**♦*»#**»***#**#♦****#**#♦#*»*#»****#***♦**********»***,# 

CONTINUE 

«' 

UVC  =RHOT  t*VDI  K.  )  *ABSI-  lVLHM)*CD!Kt)  /Kh/v  l  K  )  /  KM  K. ) 

DVIK  )  =DVC-DV ( K ) / 2 • 

ENMPI I )=ENMP{ I )+DN(Kl*DMC 

TUMP  =  SUMP  +  DTCTT7  «  ( vro-U)  *  TJX 

SUMPA  =  SUMP A  +  DVC  *EM ( K ) *DN ( K ) »DX 

C 

WW ( K  )  =  W 

75 

1200  CON TTNUE 

Q  #####*#*#***##**##*###*#*#*****####*#■*##*#*#*****#*#**#**###**#*## 

RR=RG*TG 

SDoPD  =  SCJRTF  (  RR*uAM*G  ) 

IFTU‘.5TY5D5PD>120  5  »12lo 

1 203  XT  =  0 • 

NP  =  0 

1210  AX= A ( X ) 

PH  I P  =  - ( ENMP ( 1 J+ENMPI2)  ) 

PHI =FA( 1 ) +F  A ( 2 ) 

U  =  PHI/RHOG/AX 
EMACH=U/SDSPD 
EMACh2  =  EMACH*  *2 

RATP  =  (  F  A (  1  )  / F  A  (  t  )  *LNMP  (  l  )  -t/MMP  (  1  )  )  / F  A  (  2  ) 

TF(  RAT* GT •  RATMN  •  Ai^J.T*  A  T  •  l  T  •  kA  TkiX  )  i^^4*  Illc 
1222  RATP=0. 

1774 - CONTI  .‘JUT - 

C 

A 1  =  AX  -AP ( X )  *  DX 
A2  =  AX 
A1DA2  =  A 17  A  2 

PI  =  U  *  ENMP  *  DX  /  (0*A 2) 

- P2  =  -SUMP  !  (G#A7T~ - - -  - 

P3  =  -SUMPA  /  (G*A2> 

7A  "= -"RFfOG'"1r'tflF*‘’2 "  ■X  (  A 1  /  A  2  -  1  •  Q")  /  b 
C 

—  p-  -  =  -pr  v'RI  +72'  +  D2  +  P4 
C  - 

- *  1  FT  KK ICKK-TC  )  9  89 1 ,  R  5  ?  2  ,  9’ 6  4  2 - - - - - 

9891  KKKkK=KKKkk+ 1 

'uo'tc'ws - 

9  b  92  WR.I  TE  (  6*8009  )  SUMP ♦ SbMPA . P  , P 1 , P2 . P3 . P4 .  RhuG.G.Ai ,  A^»  muAc 
8  007  '  FORMA' T  T I  h'-  76  FI  3  .  a  /  6  F  1  3  *  o  ) 

<<<.<<  =  ? 

9893  CONTINUE 
C 

X  =  X  +  DX  2~ 

L 

001260  I=T,2 

FA  I  I  )=EMDOT(  1  ) 

Ti^E'ETTT) - 

001260  J=1*N 
k  =  J  f J-2+ I 

I  F  (  EM  {  K. )  )  1260*126.  .1253 

c 

12  3  3  EM  (  k  )  =EM  (  k  )  -+QX*OM(  k  ) 

FA  I  I  )  =  F a  (  F)  -ON  (  k  )  * EM  (  \  j 
TL<  K  )  =TL( k)+OX*OT ( s ) 

V  (  K  )  =  vTkf +  DX  *  D  V  (  K  ) 

c 

I F ( t  M  I  K  }  ) 1180.1  180* lcbu 
118b  EM ( K  )  =0 • 

TL  (k)=C. 

V  (  K.  )  =0* 

RSI K  )  =  0* 

K  )  =0. 

DiNT  l  I  r=DNf  (  1  >-DN(  K  ) 

ON I  K ) =0 .  _  7$ 

i c 6u  CONT I N JE 
C 

K  A  T  =  F  A ( 1 ) /FA ( 2 ) 

IF(RAT.LT.RATMN)  1261.  ic62 


1261  RAT  «  RATMN 

. TSO“TO“"TnW . 

1262  CONTINUE 


1262 


•  U  I  • K  M  | HA  I 

1263  RAT  *  RATMX 

'I76TT - "CUflTTNUF - 

GO  TO  104 

'■TOVTfP"s-TJ - 


MU 


*******#***«****#****#**#**###*#««##«***#********«»*•«*#*«#****«*« 


RAN  =  SORTF(AX/3. 14159) 

- ■ET7r-=“7a“7’"7fXMTTT - 

IF ( RAT  -OFSTOC)  31»31»32 


Si.  f\V_vnD  -  “Cnilr  l  i  j  w  f  i  •  j 

FCOMB  =  F A ( 1 )  *  (1.  +  1. /OFSTOC)  /  LMD 
- T70“7O"T3 - 

32  RCOMB  =  -ENMP ( 2 )  *  (OFSTOC  +  1.)  /  tMD 

— rctriqB’*v-rAi7r-*"TT;vof‘5rorr'7‘TPii5 - 

33  RLCOMB  =  RCOMB  *  RAN  /  ETA 


PO  =  P*( TSTG/TG )**GAM5 

- Trcr=TTH0*vs“-ffSi5§PD_ - 

FVAP 1  =  FA ( 1 )  /EMDOT ( 1 ) 

- 'FV7CT,7*FAl"2y7EWDUTT2T“ — 

FVAP=  PHI/EMD 


=-ENMP ( 1 ) /EMDOT ( 1 ) 
RVAP2*-CNMP ( 2 ) /EMDOT ( 2 ) 

rv7Tp*~phip7emd 

DO  320  1=1*2 

savftd ) =0. 

SBVR32 ( I ) =0« 


SVR(  I  )  =0. 
SVR32 ( I )  *0* 


SR(I)*0. 

SR 3  (  I  )«0. _ _ 

N-NSET (I ) 

DO  320  J*1#N 


* J+ J-2+I 
IF(EM(K) ) 320*320*310 
‘3*  1 U  ENR*DN(  K  f*R~s7Kl 


SR  3  (I )*SR3( 1 )+ENR*RS(K)*RS(K) 


SVR( I )*SVR( I H-ENVR 

ENBVR»BETA( K)»ENVR _ 

SBVR  I I ) -SBVR < I )  -ENBVR 
R12*SQRTF(RS(K) ) 


SVR32I I )«SVR32( I )+ENVR#R12 
SBVR32  ( I )  *SBVR32  ( I  )+ENBVR*R12 


324  AB(JZ)  *  X 


0R( JZ»2)  -  FVAP2 


m: 


0R( JZ»4)  •  FVAP 


1:17.1 


0R( JZ*o)  «  RVAP1 


OR  t  JZ  *  7 )  =  RVAP 
uR ( JZ  » 6 )  =  LMAcn 
OR(JZ.49)  =  U 
UK  t  JZ  » 8  0  )  =  FCOTT5 
OR ( JZ  » 8  1  )  =  RCOMB 
OR ( JZ  » 18  )  =  P 
OR ( JZ  #  28 )  =  PO 

OR  ( JZ  *36")  =  TG  -  -  -  ' 

OR ( JZ  *48 )  =  TSTG 

J2  =  jZ  +  1 - 

328  IBZ  =  IBZ  +  1 

WR  I  T  b  (  6  ♦  3  3 0  )  <  T  1  T  L  b  ( T  )',T=  1*  397 
33o  FORMAT ( 1H 1 / ( 1H  13A6)) 

WRTfE'f  6*  3"40Tx'»'OVTG  .  TSTG'tP7PO#RAT'.'EMACH 
34U  FORMAT  (  3HOX=F6.3»4H  U=F7.1,4H  T=F4.0»4H  TO-F4.0.4H  P  =  F7.<i, 

I  4H  PO  =  F  7  •  2  *  6H  0/fr  =  F7.3»7H  MACH=F5 •  3  ) 

WR1TE(6.350)FVAP1*FVAP2.FVAP 

35o  format T 26  h  vaporized  Fraction  '~o~F9.6 *Vh'“"  fTv . t> » 7h  DOTHhVFb ) 

WRI TE (6,360) RVAP 1 .RVAP2.RVAP 

360  FORMAfl^BH"  VAPORIZATION'  RATE'/ I N  “  '0F9~.6.4H  FF9  • 6~»7h  ""cothbV.b  ) 


N=NSET  (  1  ) 


WRITE(6.370) 

3 7o  FORMAT (70H0OXIDIZER  DROPS 

R(MIL)  MASS  (  Lb )  VUN/StC)  U-V/A  TbM 

IP  NUMBER  ) 

DO  380  J=  1  »  N 

<= J+ J-l 

I  F ( EM ( K )  J  371 ,371 , 375 

3Tl  RPO  =  0.0 


VDOUT(K)  =  __0.0 

""  CaO  'f C)  377"  . 

375RPD  =_DM(<)/EMINT(<) 

3  7  7  "  F'R  A'C  ~T~  "  '="  ~  "( ~E  mT  N~f  I  )  ~  ~  -b  M(7  f )  7  E  M  I~N  tT  Rf 
MNM  =  J  +  8 

OR(JZ-l.MNM)  =  RS ( <  1  *  l.E  +  3 

MZM  =  J  +  28  _ 

0R-(  JZ'-T.'m ZM)"'=  ab s'f  (  v'do'ut 'Tm 

MZN  =  J  +  49 
0R'(’jZ-Y I M Z N  r ’="  V (<T ‘ 

3  80  WRITE(6,39Q)J,RS(K),EM(K),V(K)  »VDOUT(K,)  ,  T  L  (  K  )  ,  DN  (  K. )  , FR ACT  ,  KPD 

390  FORMAT( 13X. I  2 » 3PF9 . 4 , OPE  1 3 . 5  »  F5t0,F9.4,Fd.2,Lli.3»4X, 

1  Ell.  3. 4X. Ell. 3) 

WRI TE (  6.400 )  ' 

400  FORMAT ( 12H0MJEL  DROPS  ) 

""N=NSET(T> 

_  DO  410  J=  1  >  N _ _ 

K  =  J  +  J 

I  F  (  EM  (  K  )  ) ':01.401,40b 
'40T"R"PD"  ="  cuo 

VDOUT(K)  =  0.0 
'GO  TO  *  407 

405  RPO  =  DM( K ) /EMI  NT  (  K. ) 

ZfTJTT-FRACT  =  '( EWTTTTR )  ^EM<K))  /EMI KlY  <  K ) - 

MNM  =  J+18 

ur  rjz’-y.MNfo'r  ='  Rs<  '#T.i+r" 

MZM  =  J  +  36 

URT JZ^T » M'Z W r "=  AB'SFTvUOUT(XTT - 

MZN  =  J  +  59  7, 

- ukuz-i.mznj  =  vin -  - 

410  WRI  TE  (  6.391  )  J.RSdO  .EMU)  .VU)  .  VDOUT  (  K. )  »  T  L  (  <  )  .DNU)  .FRACT.RPo. 
lXDTSUTX-y - 

391  FORMAT ( 13X. I2.3PF9.4.0PE13.5  ,  F5.0 .F 9.4.F8. 2 • t 1 1. 3 »4X . 


1  Ell. 3. 4X, El 1. 3.  1  13 ) 

UO'‘4TO*  1-1*2  . . . . — . .  "  . . 

K 1 C  (  I  )=$R<  1  )/DNT(  I  ) 

- KTDT'l  )  =  1  e  R'j  11  )  /  l)N  1  I'l  ))  »♦  1  HKU -  - 

RV321  {  I  )  =  (SVR32  (  I  )  / S VR  (  1  )  )**2 

R6\r3?irry=i5nvR  J2  (rr7T6VRTm*"*2:"  — . "  ' 

4Ju  CONTINUE 

'  ‘"'rALL'TERP’  1 K  1 0  *  RS  ♦  vUOUTVNlTET^VDrO') - 

CALL  TERP  ( K30  »  i-i S  ♦  VDuLJT  »  NSt  1  *  VD30  ) 

CALL  TEkP  <  RV32 1 »  Kb  •  vuuu I  ♦  NaE T , VD3  2 1  ) 

CALL  TEkP  IkcV321  »Kb»VL>OuT  »NSt  T  »VL>d321  ) 

RTIC  "='TsKn  )  +  "S'KTYTT“/r_TuNT*ri'r"+" BnTT2T) 

KT  3  U  =  (-(SKi(l)  +  SR  3  (  2  )  )  /  (DNT(l)  +  oNT  (  2  )  )  )  **ThKL- 

'  ■TfTV3rr'="TRV37rnT'+'TfmTrrn"7r"2v - 

k  T  d  V  3  2  1  =  (  RbV  32 1  (  1  )  +  RHV32K2))  /  2. 

- VZTTZ  =  rCToTT  rPHTTOTH + — VD3U 1 2  )  *R 30  (  2  77  /  ( KTOTT)  "+  k  j  0  CZ7 ) ~  — 

VOTlw  =  ( VD10 (  i ) *R1C ( 1 )  +  VD10( 2 ) *R10( 2 )  )  /  ( K 1 0 (  i  )  +  kiO(c)) 

"V3T3rr*=  rvuniTn  *7732 it  it-  + '  W3  2 n7r*W3  21 17 ) >  7'tk  v32i  ( 1 »  + 

1  k V 32 1  (  2  )  ) 

- V0T572T"  V-  CV  76  37IT1T*TT6  V37ITTT '  VU572TT' "2T*R  BV721T2TT ' ~T 

1  (  R tJ  V 3 2 1  (  1  )  +  RdV 32  1  (  2  )  ) 

lJoh7TJT=TT2 

REN (  I  )  =  RIO ( I )  *  RcC 

TTENTT+T r'^TfoUV ’*  RK - - - 

REN( 1+4)  =  RV321  (  I  )  *  REC 

- 7E7^TT+¥r'"="RBV371TIT' * "FTEr - -  - - 

403  CONTINUE 

KtiM  (  9  )  =  HT  10  *  ktC 
REN (  10)  =  RT30  *  REC 

R£N(TrrV'Rfv32T"*'RtC 

RENI  12)  =  RTBV321  *  REC 

---- ■-  ^ .  . 

3D2  =  8HFUEL 

dOl  =  8HOXICIZER 

dRN(l)  =  RVAP1  *  RAN  /  ETA 

d'RN<2  )” ="” RVAP2  ”*"  KAN~7‘t“f  A-  . 

dRN ( 3 )  =  RVAP  *  RAN  /  ETA 

-mrjz-ue2~rv  ritnytt 

OR( JZ-1»83 )  =  REN ( 8  ) 

0R(JZ-1,84)  =  RtN ( 1 2  ) 

OR ( JZ-1 ♦ 8  5 )  =  RLCOMb 
OR"(  U  Z~- 1' .  8  6  )*"="  8*RN73  ) 

OR ( JZ-1 »  c  7 )  =  AdSF (VUb321(2) ) 

■cRTjr-T»FbT  "='_A'ds>'(VDb-57rrrrr 

CR ( JZ-1 » 89  )  =  AdSF ( VOT8321 ) 

- I F  ( JZ- 1 3  o ')  7  j  7  5 »  7JT5TT59 - 

7373  mR I T  E ( 6  *420 ) 

7"RTT1""To'»”4  3"c Td D’l  *T<  l  c  ( YT Vv dToTIT *  R L n  I "l”)  Vr  3'o"(T”)"»  vd  30  ( T )  *  k  t  n  (  3 )  * 

1  RV321  (  1 )  *VD321  (  1  )  *Rc.N(  p  )  .R8V32K1)  » VDd32 1  ( 1 )  *KtN (  7  ) 
w'R'l  T  ETs".  4"30"rBD2VR70T27^v'0rcl'2'r*“RE'N72T»R3b72T*’VD’3'6  (~2  )  *RfcN("4  f  * 

1  RV321  (  2  )  t  VD32  1  (2  )  »KbNl6  )  *RbV321(2 )  ,VDd32i<  2  )  » KtN ( 8  ) 

- WRITE  (6,430  oD3.«l  1 0  7  VUT  lOTITEKi  I  9 J .  RT30  .VDTO.kENT  10 )  ,kT\/32  1 , - 

1  VDT  321 ,REN ( 11) ,RTbV32l » VDTd32 1 »REN ( 121 
“4  2  u  "TOTTm  a  'TTTho*  137  #  6hT?  T  m  I  IT7Tx  » 4HTJE  l"v“8" x73  hREd  1 
430  FORMAT! lH0,A6/( 14X , 3PFb.4 , OPF 12.4,0PF12.4/ )  ) 
w RTTYT 6" , 4 5‘CTTb Rn ITT  *  T=TT5 ) 

44u  FORMAT! 1H0*13X,23HBURNING  RATE  PARAMETERS  / 14X * 8H0X I DI ZER  , 

- 1  F1I\2.4hFUEl — ,  F11.2,5MTCT4L  /FTT'.fj - 

WRITE  (6,441)  RC0M3,RLCoMd,FC0MB 
44 1  TQRR7TT  TIhu ,  I6h7JmdUS TTuTj  RaTE  rFTrSTTXTTSH^UMF^  PARAMtftK 
1  5X , 14HFRACTI0N  COMb.  ,F9.6) 


79 


WR I TE ( 6  » 55 )DX»GAM 

55  FORMAT  (  18H0INTEGRAT10N  STEP=F7.4,5H  I  NCH5X6HO  AMNIA  =  Fi>  •  J  ) 
GO  TO  250 

'999  dCDY  1 1 )  "s’-ffHGTF  RAT  I - 

JZ  =  JZ-1 

KZ  =  AB( JZ )  +  2.0  ' 

UPA  =  KZ 
XUU  =  UPA' 

UPAM  =  UPA  -  1.0 
DO  990  I  =  l.JZ 
K  =  I 

I  F  (  AB  ('<')'  -  T.O  )  990  *  99'oV  99  1  "* 

990  CONTINUE 

^  ........... 

991  M  =  K  -  1 

- DO  993  T"=  T7R - 

993  AB( I )  =  AB  (  I  1*3.5  +  1-0 

frrjr  -----  996 ,"996"'~994 

994  DO  995  I  =  K  ,  JZ 

'"AFTTT""='"r  ABITI-IVOT?  IJPA'R*  3V5"+‘ 4;“5 - 

996  BCDY ( 2 )  =  8H0  OF  CUM 

- BTTTrTTJ—^HHFr^TrSTS - 

6CDX ( 1 )  =  8  HD  I  ST  ANCt 

- srDxrrr'="bH"iiT'TNrH - 

dC  DX ( 3 )  =  6HES 

- TD"=-T*7rrwx - 

UO  =  ID  +  1 

call  gkaPH  (  Ab«G,sx ( l .  l )  »xuu,uo  ,  BCD* ,  BCdy  ,  jz  » l ) 
dCDY(l)  =  8HFRACT  1  ON 

- oCUYT2T'="Bh- vATJ0T<TZ - - -  ' 

BCDY  (  3  )  =  8HED 

- -3r-(5yr4-r-=--1 m -pxf£L" - - - 

BCDX ( 5 )  =  8HOX I D I ZtR 
dCDX ( 6 )  =  8H  TpTAL 
BCD_X  (  7  8_HC_OM6U_ST__ 

CALL  GRAPH(Ab»  OR ( 1 » 2  I  * -XUU » - 1  •  0  »  dCDX,  dCut  *  jz,  j.) 
CALL  G_RA_PH_(_Ab.  OK  (1 . 3  )  » -XUU  » - 1 . 0  »  oCDX.  dCDY,  JZ.-Z) 

C  A  L  L  GRAPrll  Ad  »OR ( 1 »  4)*  —  XjU,-l»G»dCDX*DCDf*  j  z  *  ~  i  ) 

CALL.  GR APh ( AB » OR ( 1 » 6 0 ) » + XUU ,  —  1 • 0 » dCDX * dCD Y »  JZ»  -4) 

- 3rDV(1)  =  BHvAFCR  1 ZA - 

act:  Y  (  2  )  =  8HTION  RAT 
BCCY'(  3')"'="8'HtT'  ' 

UO  =  0.0 

- DO-W"r=''lV"JZ  . 

831  UO  =  MAXI ( JO  »,  OR (  I  »  5  >  »  OR ( I  *  6 ) .  OR  t I  *  7 )  ) 

UO  =  UO  +  1.0 

CALL  GRAPHIAB,  OR ( 1 *5 ) * -XUU.-UO  »  dCDX  »  dCDY,  JZ  »  + 1 ) 

'  CALL'GR a'phTa'b 7"or"(  liV) V-xuu."-u6'  »"bcdx',  BCDY .  JZ » -  2  ) 
CALL  GRAPH ( AB  »OR ( 1 »  7),-XUU,-UO  ♦BCDX.dCDV.  JZ  »  -3) 

CALL  'GRAPH  (AB  *6r"(T  *  8_I  )  »  +  X  UU  V-  UO  ”  V  BC 5  x" *  BC  D  Y  .  JZ  .  -4  ) 

BCDY ( 1 ) ,=  8HGAS  MACH 
BCDY  (2)'  '«=  *BH  NUMBER 
dCDY  (  3  )  =:-'8H_ 

CALL'GRAPH  f  Ab  ."OR  (  lTsT,  XUUVTroV  BCD  X  ,  dCDY,  J  z*  1) 
dCDY(l)  =  8HOXIDIZLR 
"3'CdYT2T"=''8H  .drop'ra" 

BCDY  (  3  )  =  6Hfc  (MILS) 

UO  =  0.C*  V. 

KZ  =  NStT(l)  +8  80 

'DO'  '  TO  O'C'T  '  "=  '  T  J  Z 
100U  UO  =  MAX1F ( UO»OR ( I *KZ ) ) 


I 


00  1015  I  =  1,JZ 

- I-j1-5-{jR-rr--9Irr"  "0RT1V4T1 - - - - - 

MZN=U0 

- UO  =  MZM-T - 

NMN  =  NSET(l)  -  2 

rAUL' oR AP'r:i  asvo"< tt;  ~9  >  r-xiiu, uo  ~.&'cdxVbcT)y  ,"jr . + i  i 
DO  1020  1=1.  NMN 

- - 

102o  CALL  GRAPHt  Ab»OR( 1 »NZ ) »-XUU»UO  , BCDX . BCDY , JZ . -2 ) 

- CALL  GRAPH  (  Ad *GR  ( 1  »KZ  )  »+XUU»"GO — » bCDx » BCDY , JZ  » -3  ) - 

BCDY (  1 )  =  6h  FULL  DR 

. . "H‘CI5YTZT-=’”B'H0'P"TrA'LJ'T - 

BCDY (3)  =  8HMILS) 

- - 

KZ  =  NSE  T ( 2 )  +  lo 

- DO  TU7D  I  -  TTJ7 - 

lu5u  UO  =  MAX  IF ( UO  *0R ( I » KZ ) ) 

- ffZPUD - 

UO  =  MZN+1 

- NWj_.._-N^1Err7.r____7. - 

CALL  GRAPH(AB*OR( 1,19) .-XUU.UO  , BCDX , BCDY , JZ . +  1 ) 

- DO'T070 - 1  -  1 «  NMN - 

NZ  =  I  +  19 

- ITTrj‘rALT""GTrAPHi7vF.“crRTT;TMZT?-‘7U0",i;;cr‘“TBUD'X‘;‘FajYrjr?-‘zi - 

CALL  GRAPH( AB»0R( 1 ,KZ ) ,+XUU.UO  ,  BCDX , bCDY , JZ . -3 ) 

. . — NW3--=“Nsrrrrr - 

bCDY(l)  =  8H0XIDIZER 

- baJ77 zy  —  Hh '"RE LI'  MA - 

BCDY ( 3  )  =  8hCH  NO. 

. 

KZ  =  NSET ( 1 )  +  28 
_____________  _ 

_  DO  1 C75  I  =  1,  JZ _ 

1U75  UO  =  MAXI  (UO,  100. *0R ( 1 , KZ ) ) 

UO  =  <U0___+_  2_.  0_)_/_l  00  • _ 

c'al_l”gr_aph ( A6 .or  ( i", "2 9  r » ”xuu7uo • bc d x 7b~CDy7  J Z  7T) 

DO  1080  I  =  1,  NM1 

_ 1080  CALL  GRAPHtAB.  OR ( 1 , NZ ). -XUU  .  UO , BCDX . bCDY , JZ .  -2) 

CALL  GftAPN ( Ab ♦  OR(l.KZ).  XUU » UO » bCUX , bCDY » JZ »  -J ) 

_ BCDY(l)  =  8H  FUtL 

NM2  =  NS"ET(2T'-2  ' 

KZ2  =  NSET ( 2 )  +  38 _ 

U0___=_0.0 

DO  1085  I  =  1 »  JZ 

Tu53  CJG  ^  maxkuo,  igo.*or ( i ,kz2 i j 
_U0  =  (UO  +  2.01/100. 

^Qr"^fi7rp»TrAii7_oi?(Tr3r)_;7xi7u*_u67"‘B£'Dx*"bcrDY7j_z7'rr 

DO  1086  I  =  1.NM2 

. 

1086  CALL  GRAPHt  AB.ORQ.NZ)  »-XUU»UO»BCDX, BCDY, JZ. -2) 

- acrcwm  ab  ,or  i  1 ,  kz2  t .  kuu  tout  BtcxTUCD'Y  ,-jz-,  -3 > - 

KZ  =  NSET ( 1 )  +  49 

- ■5CT7?TIT“=__l!T!CrxlUlTETr - 

BCDY ( 2  1  =  8H  DROP  VE 

- -5Cff?T3T__=__OTU_TN?-5Er - 

BCDX (4)  =  8H  GAS  VEL 

- DO  677  1  =  1,9 - - - - - - 

KK  =  I  +  49  at 

- -  - 

K  =  J 


— J 


I F ( OR { J.KK)  )  676.  677,  676 

o7o  continue: 

K  -  K  +  l 

677 - LOB ( I )  =  K  -  1 - 

UO  =  0.0 

DO  1088  I  =  1.J2  - - 

UO  =  MAXlF(UO.OR( I *49) ) 

nzz  =  uo/rooo;  +  1.0 - - -  " 

uo  =  NZZ  *  1000 

M'NM'  =  NSh  I  (  1  )  ~  I - 

CALL  GR APH ( Ab  » OR ( 1 » 49 ) » -XUU  »-UO *  BCDX  »  oCD  Y  »  JZ  »  1) 
DO"’T"OF"9  --1-s-ivmNR- 
NZ  =  I  +  49 

KiT'^'TOBm - 

1089  CALL  GRAPH( AS.OR ( 1 .NZ ) .-XUU.UO.BCDX , BCDY .KK.  -2) 

- fX-r-nJFTTJ - 

CALL  GRAPH(AB*OR{  1 » KZ  )  » XUU .UO.BCDX .BCDY » KK»  -3) 
-Sr0yryr-£--ffFr— jtqtt-q- - 

DO  68?  I  =  1.  9 

- 

DO  686  J  -  l.JZ 

- - 

I  F ( OR { J.KK)  )  686.  687.  686 

■£B7r"ctrfiTrNor - 

K  =  K  +  1 

"'758T  U0BTT7""-'X"“-”T - 

NMN  =  NSET ( 2 )  -1 

- KZ  '  TTSET T2  J  ‘ "+ — 5V - 

CALL  GRAPH(AB.OR( 1.90) .-XUU. -UO.BCDX.  BCDY.JZ.  1) 

. . "0U"T09TT"=— f.NMlj 

NZ  =  I  +  *9 

. . .  <X'="CObTn - 

1091  CALL  GR APH ( AB  »  OR  (l.NZ).-XUU.  UO  »  BCDX  »  BCDY  »  KK » -2 ) 

- =  L0b(9> - 

CALL  GRAPH ( AB .OR ( 1 » KZ  ) ♦ XUU .UO .  BCDX.BCDY »KK.  -3) 

BCDY  (  2  >  =_  8hYN0LD_S_N 
d COY  (  3 T  =  8H UM BE R 

BCDX ( 4 )  =  8H0X I D I ZtR 

bcdx ( 5 )  =  8h  Fuel 

UO  =  0.0 

-jo”  j  y  I  o-  T  ‘  T*“  j“z 

1110  UO  =  MAXI (UO»OR( 1 .82 ) .OR ( I .83) »0R( 1 .84) ) 

- rr[UO-”-Yoo5'.'r"6‘5’6V‘656r657 - 

656  UO  =  1000. 

GO  TO  658 

657  NZZ  =  UO/IOOO. 

. 

658  CALL  GRAPH  ( AB .OR ( 1 » 8 2 ) .-XUU . -UO . BCDX . BCDY , JZ .  1) 
rALT"5ITA¥>H‘TABV0ftTIV8y)V-xuu%"W,'bCDxrbC0Y",j27-2)' 
CALL  GRAPH  ( AB .OR ( 1 ♦ 84 ) . +XUU .-UO. BCDX » BCDY » JZ »-3 ) 

- BCDYTTT  =  'BRBCIRNITIG - 

BCDY < 2)  =  6HRATE  PAP 
■dCT5YT5T‘=  'EThAmETE'R' 

BCDX ( 4 )  =  8HC0MBUST 

BrDxiTr'=‘^FiTOTAr - 

UO  =  0.0 

- 7TO  TT7T  I  =  TTIZ - H - 

1121  UO  =  MAXI  (UO»10.*OR< I .85) .  10. *OR ( 1 . 86  )  ) 

- IFiUo"-"I7r;r"B'?iv'‘S’2yr"s^7 - r - 

821  UO  =  (UO  +  1.0J/10. 


^521 
c  24 


GO  TO  824 

■— cro"=‘'ruo'+_TDvr/iu; - - - - 

CALL  GRAPH ( Ab *  OR ( 1*83)  •  -XUU»  -UO.  dCDX.  dCDY.  J 2*  i) 


'-MLL  U  r\  Mr*  n  l  Mu  »  UK  l  i  *  0  fc  J  »  +  A  UU  .  ~UU  »  DO  l>  X  *  D  O  U  T  *  J  L  »  ~C  ) 

dCDY (  1  )  =  onSTAd I L I T 

‘FCDTT2T"=‘"5hY”DT?0'P"U - 

bCOY ( 3 )  =  8HELTA  V 

'FCUXT4T~=“~8TTT[jEL - - - - - - - 

dCDX ( 5 )  =  8HOXIDIZER 


Hit] 


"IT3T 


DO  1131  I  =  1,JZ 

”--u'c”=_'WA‘Xi“TTioTiTiov*uRTrT8TTviu‘0'.TFoi?TT;YaT;TO~o';:«f'OKrr»"B'9rr 

00  =  (00  +  2.01/100. 

■OLT"GR’APVrTA'B",URTIT8TTV-yOU;-TiO'._B"CDy;'ffaJ?TOZr'IT - 

CALL  GRAPH  ( A6.0R ( 1 .88 ) .-XUU.-UO.BCDX .bCDY. 02 .-2 1 


oCDY(l)  =  8HCHAMBER 

•FCUYT2T"V"5HFT7E^5URr' 
8CDY ( 3 )  =  8H 
••3-CUXT4T“=‘-FK5TAT"FRV 
dCDX ( 5 )  =  8nT0T .  PR. 


ititi 


DO  801  I  =  1,JZ 

"80T'TJ0"'='"MA"xrnjc7c5RTT"ri'8')V0Trn"r2"5rr - 

NZZ  =  UO/IOO.  +  1.0 

- U0"="fi7Z*-rTn; - 

CALL  GRAPH ( AB »  OR (1*18)  *  -XUU,  -UO.  dCDX  »  dCDY »  JZ »  1) 


H(Ad*OR(  1.28)  »XUU»-UO»dCDX»BCUY* JZ.-2) 
dCDY ( 2 )  =  8HTEMPERAT 


¥cdyT3T'="-8hure' 

dCDX  (4)  =  8HSTAT_T. 
dCDxT5T"="8HTCTV 
UO  =  0.0 


0  802  I  =  l.JZ 
802  UO  =  MAXI (UO,OK( I .38) ,UR( I .48) ) 


uzz  -v-uo/'icroo  .""vt:  0 

UO  =  NZZ* 1000 

TA'LL"GTrA¥>"H‘TABVORTI  ,"r87V-Wr-U07BCDxrd"CDY7jzri"")' 
CALL  GRAPH  ( AB .OR ( 1 , 48 ) >  XUU.-UO.BCDX.BCDY . JZ »-2 ) 


SUBROUTINE  GRAPH  (  ABS.ORD  *UPA  »UPO»bCDX  » BCDY  *NP »  NO 
DIMENSION' AB'SUOC)  » ~QRD(  150  F»  BCDX  (  7  )  »  BCDY  (  3  )  *0(400) 
DATA  (JZ*G) 


IF(JZ)  9.8,9 
CALL 'PLOfSTG *225  » 1 7 F 
JZ  =  1 

’  ITC  AX  tf  ”T  N’C  OWN  G“  “  U  AT  A" 
NPT  =  NP 


*  XABSF(NC) 
APA  =  ABSF(UPA) 


UUlH  *  1  (-V 


•  Af>0'"*— ABSFTUPO)  . 

DO  14  I  =  1 »NPT 

0^07rr"='‘ABsrroFDTIT7UP09.‘0>"+"o75' 

IF(NC)  7 1 »  3,  3 


UR  ORDINATE  LAB 
I F ( APO  -  10. )  124,  124,  22 


LAYOUT  ORD  LAtJtLS 


W'=T 

DO  125  I  =  1,3 


i“T-"r - 

I F ( APO  -  1.0/10. 0**N )  125,  125,9126 


Ut 

YB  =  1 . 0/ 1 0.0*#N 

- 

I F ( APO  -  3.0/10, **N )  134,  134,  135 

---YB^'TBVZVtr - 

YZ  =  9.C/AP0*Y6 


.0/ 

GO  TO  25 


9126 


"MT'="2 

IF (APO  -  20.  )  24,  24,  23 


■rrfAPO'-TOOOO  126T'26',"""2T"' 

YB  =  100. 

1 F ( APO  -  300.)  142,  142,  143 
YB  =  YB/2.C 


YZ  =  9.0/APO*Y3 
GO  TO  25 
'Y6""="2C0. 

YZ  =  1.8 


GO  TO  25 

I F ( APO  -  2000.)  28,  28,  25 


YB  =  500. 

YZ  =  2.25 
"'GO''TO"25 

YB  =  1000. 

YZ  =  9. 0/ APO*l 000 . 
GO  TO  25 


YZ  =  9.0/APO 
Yd  =  1.0 


XLA  =  0.0 
XAB  =  1.0 
X  =  XAB  -  .07 

CALL  NUMBER(X»0.4,.07»XLA,G.0»2) 


CALL  PLOT ( XAB »  0.5,3) 
CALL  PLOT (  X  AB »  9.5,2) 

— xA¥"V-xaB"V'.8?1>"" . 

CALL  PLOT ( XAB ,  9.5,3) 


LL  PLOT ( XAB ,  0.5,2) 

XLA  =  XLA  +  .25  14 


X  =  XAB  -  .07 

CALL  NUMBER(X,0.4,.07»XLA,0.0,2) 


XLA  =  XLA  +  .25 

- JTAF'V-y-Ao  -V-VBT5 - - 

X  =  XAb  -  .07 

l XLA  -  1 • 0  >  41 1  4 2.  42 
42  CALL  SYMo0L(3.0. 0.2,0. 14, 6CDX.  0.0*  24) 

TATU"NUMdERiyvjr^;voT;xL7r;F:ir,'2‘r - 

CALL  PLOT ( XAB ♦  0.5.3) 

FAL“C~"PL'0TTX'AB~."5i5»'77 - 

I F ( APA  -  1.0)  41.  41.  45 

55 - Z  =  3.5/MAPA-1.0) - 

XAB  =  4.5 

TFT7'~-“VBT5T~5&T~57T~57 - 

46  Z  =  2 » 0#Z 

— yi7TVT;o' - 

XB  =  2.0 

- GOTUT8 - 

47  XLA  =  2.0 

- ¥'="T.T - 

48  XAB  =  Z  +  XAB 

- rF[xAD""-""S'.'0'O0o'IT'4?T'497'r4 - 

49  X  =  XAB  -  .07 

50 — CALL  FL  <3'T  T  KT6T"  '9TF7Tj - 

CALL  PLOT (XAB.  0.5,2) 
•CATL'NUMB^"x7o":4";76TrxTA7u."o,'27 

XAB  =  XAB  +  Z 

- rFrxAF‘-"B".*o_(yo5rr5r,"5r,"55 - 

51  X  =  XAB  -  .07 

XLA  =  XLA  +  XB 

52  CALL  NUMBER ( X .  0.4, .07.  XLA.  0.0.2) 

CALL  PLOT (  XAa,9.5,  2) 

---- — yg  . 

GO  TO  48 _ 

54  XAB  =  XAB  -  Z/2.0 
_ X  =_XAB  - _• 07 

xITa  V“xl~A  "-"xB/2’."o“ 

I F (  XAB  -  8.0)  50,  50*  58 

55  XAT"="'xa’§""-'z/2".0 

_ X  =  XAB  ~  .07 _ 

XLA  =  XLA  -  XB/2.0 
_ I  F_  (  XAB  -  8.0)  52*  52,  59 

56  "7z"V"-YZ 

Y  LA  =  APO  _ 

vd”V""-9a""+‘  76oo‘5T5  . . 

_ VAB  =  9.5 _ ; _ 

AZ  =  0.0 
A  «  -i.O 
GO  TO  60~ 

59 _ YLA  »  0.0 _ 

Y~A6  =  0.5  •  . 

AZ  =  9.5 


6u  CALL  PLOT  (  p.O.YAB*  3) 

- -CTnrrajrr-  rTcrjYABTTr - 

Y  =  YAB  -  .035 

- tftttzt  roT»~~rcrzrr-ro'5 - - - 

103  call  number <  o « 7 »  y  *  .07,  yla.  0.0.  3) 

— — GO  TO  106  - - - - - 

104  CALL  NUMBER ( 0. 7  » Y  ,  .07.  YLA,  0.0,  2) 

- GF"T0”I55 - 

105  I F { MZ  -  1)  108,  108,  109  »$ 


10b 

109 

CALL  NUMBER (0.7. Y.  .07.  YLA.  0.0.  1) 

GO  TO  106 

CALL  NUMBER (0.7. Y,  .07.  YLA.  0.0.  0) 

YAB  -  YAo  +  YZ - 

IF ( (AZ  -  YAo) /A)  70.  62.  62 

62 

Y  L  A  =  YTA  +  Yo 

Y  =  YAd  -  .035 

IFfWZl  113".  H4.  1T5 

113 

CALL  NUMBER ( 0.7 » Y.  .07.  YLA.  0.0.  3) 

GO  To  116 

114 

CALL  NUMBER ( 0.7.Y .  .07.  YLA.  0.0.  2) 

- 

GU  TO  116 

115 

IF(MZ  -  1)  118.  116.  119 

HIT'  ~C7CEU~NUMBE717V7»TV~.1T7T-YC'A‘»"~0T0V"T)‘ 


GO  TO  116 


TT9 — CALL  NUMBER  ( 0  •  7  *  Y  *  .0  T~i  YLA  *  0.0.  0) 


116 

call  plot ( i .o »  yab.  3) 

CALL  PLOT!  8.0.  TAB.  2) 

YAB  =  YAB  +  YZ 

TLA  =  TLA  +  T3  ~'r" 

IF( (AZ  -  YAB)  /  A)  70.  60.  60 

To 

C 

DRAW  CURVES  ON  GRAPh  WHERE  NC  REPRESENTS  NUMBER  OF  CURVES  ON 

'"C" 

TnlS  PLOT 

71 

call  PLOT ( ABS ( 1 ) »  ORD ( 1 ) . 3 ) 

r F  f NCT-l )  72.63.72 

b3 

uO  65  J  =  l.NPT 

! 

GO  TO  91 

”77* 

— rrrN'cr"-'2T'75_T‘7Tr'7B_ 

73 

Do  74  J  =  l.NPT 

. . 7 V  "“CAT  IT  S’?  mSOI  TAB'S  ( JTV"  BRUT j~r.Vtf4  V2 V 076V-21 


GO  TO  91 

TS  IF(NCT  -  3)  79,  76.  >9 


76 

DO  77  J  =  l.NPT 

77' 

CALL  SYMBOL ( ABS ( J ) ♦ 
GO  TO  91 

ORD ( J ) *.04, 3, 0.0, -2 ) 

79 

80 

I F ( NCT  -  4)  82,  80. 
DO  81  J  =  l.NPT 

82 

81 

CALL  SYMBOL ( ABS ( J ) , 
GO  TO  91 

ORD ( J ) ,.04. 4, 0.0, -2) 

82 

I F ( NCT  -  5)  86.  83, 

86 

03 

DO  84  J  =  l.NPT 

- BTj:— CA-LU"5YM't30UrAB'ST"J)VURUTUT7.'Q"4'.'b"roror-2T 


,  GO  TO  91 

*3“  JO  87  J  =  l.NPT 


87  CALL  SYMBOL ( ABS( J ) ,  OKD(J) 

,  . 04 .6.0.0. ”2 ) 

9T  IFTUPO)  911.90,  90 

911  IF (NCT-  2)  912,  913.  914 

717  ca L L  symbol  (6.0  ,  r.5o»  .04 . 

1,  0.0,  -1) 

CALL  SYMBOL (6.3,  1.50.  .07. 

BCDX ( 4 ) ,  0.0, 

8) 

GO  TO  90 

913  CALL  SYMBOL (6.0,  1.25..04, 

2,  0.0,  -1) 

CALL  SYMBOL (6.3.  1*23,  .07, 

SCtfX  C5 ) ,  0.0, 

8) 

GO  TO  90 

714  rFTNTTT  -  31  9T5.  717,  916 
915  CALL  SYMBOL (6.0.  1.00..04, 

3.  0.0.  -1) 

— — III  II  II  !■ 

GO  TO  90 

916  CALL  SYMBOL (6.0. .T5 »  .04.  4,  0.0,  -I) 

CALL  SYMBOL (  6.3.  0.75.  .07.  BCDX (7) »  0.0.  8) 


SUBROUT INfc  COMP  ( OFR  *P  *  Y ) 

■TOMPU-STTrON"aF‘rO>TBI)STIUR"G7rSirS"_Ar'AT5r^b7fTTC"TUAM'£'Tt^PYR7rTuRt' 

AS  A  FUNCTION  OF  O/F  RATIO*  PRESSURE 


»YCOT  i  ( 6 ) »  Y  CQ2T 1(6) *YH2T1 (6) » YH20T 1(6) »Y02H (6) » 
1YN2TK6).  Y  (  7  )  » YNOT 1 ( 6 ) » YCU2 T 2 < 6 )  *  YCOT  2  (  6  )  » YH2  T2  (  6  )  *  YH20T2  (  6  )  * 
"2Y027 27 6777727 2 1  6]  »~YWr2T67 » YC0BT6T* Y'c02bT6T*  YH2 bTb")  •  YH20bT6T»"' 
iY02b_(_6)  »Y_N2_b(6UYNOb(_6)  _ 
lTAT A  TFT=  0  •  6 8"T » u  .72  5 »  f •T75  i"*"l  •  ~SZ  b  »"1.7"73»iQ.)» 

1(  YCOT  1=  U. 052,0. 043,  •  u43  *  .041*  .039».039)» 


I  Y  L 

3  (  Yn2  T 1  =  .331*  .213*  .147*  .073*  .019*. 019), 

"5T?Ff2UTY=  72B77777Yr772Tr747b67'“74"557743"0  )7 

5 (  Y02T 1=  .000*  .001,  .020,  .049,  .056,  .036), 

■6T"Y^7TT='"“.'5Bir.”;39Tr";wr”A_03",""V39"7Vr39TTY - 


776B77J 

'.72771 

•7787*7  t~S2b  ,1 . 773 

»  i  0  •  )  » 

0.052,0 

•  043, 

•  u43 » 

•  041, 

.039 

,.039) 

, 

.032, 

•  031, 

.  049 , 

tma 

•  u4p 

,  .043 ) 

9 

.331,  . 

213,  • 

14  7 ,  . 

073,  ■ 

.019, 

.019 ) , 

.237, 

.  3  7o , 

•  423, 

•  466, 

•  455 

,.430) 

9 

.000, 

.001 , 

.020, 

.049, 

.056 

,.056) 

9 

■;wr";Y03",”Y3“97 


7 (  Y NOT  1  = 

.001  , 

.006, 

.021, 

.077,  .080, .080), 

UA  1  A 

1(  YCOT  2  = 

.  036  , 

•  043 , 

•  042  , 

.040,  .033,  .033),  j 

2T?r07T7___ 

.  (J577 

.033,  .049, 

•DAY,  .542*. 042), 

3  (  YH2T2  = 

.330, 

.211, 

.  149  , 

.080,  .026»«026), 

2.TTH.2UT7___. 

.286, 

77757 

.413, 

•  447,  .44?  ,  .4777  , 

5 (  Y02T2= 

.000, 

.004, 

.019, 

.030,  .050, .050),  | 

rjsyy 


mu 


Ll«l 


7 (  YNOT2=  .001,  .009,  .027,  .083,  .019,. 019), 

- 7c-V"iT---^tro;T7„7uo7 - 

DO  10  1=1,6 

- Y.cu0Tn___r-YroTITTT.„.7r_T-1TUDT71T7„_YCUTrrTn - 

YC02B ( I )  =  YC02T 1(1)  +  A  *  ( YC02T2 ( I ) -  YCG2T1 ( I ) ) 


(I)  =  YH2TKI)  +  A  *  {  YH2  T  2(1)  -  YH2T1  <  I  )  ) 
YH20B  (  I  )  =  YH20TKI)  +  A  *  (Yh20T2(I)  -  YH20TKI)) 

Yo2BTrr"=™Y02TT(Tr"+"A  •i'YYorfrrrr'-'YoTfYrnT 

YN2BII)  =  YN2TKI)  +  A*  (  YN2T2  (  1  )  -  YN2T1  (  I  )  ) 

YnobTH  '="7 notT ( i"T"-+'‘A*(YNof2TlT-YN0TT(T)T 

10  CONTINUE 


ALL  LNGRIN  (OFR.YCO,  R  ,  YCOb,DYDX,6) 
CALL  LNGR1N  _  (OFR  ,  YC02  »R  ,  YC02B,DYDX  ,6  ) 
"'"call  lngrTY  "ToYr7yhTYr""7yh’2b"”,dydx,"6") 

CALL  LNGRIN  _  _<0FR»YH20,  R  »  YH20B  *  DY  DX_*  6  )  _ 
CTl7""ln"GRTn‘  "  (OFR  ,~Y0 2  ,'~R""7Y02B"»  "  D Y DX  ,  6  )" 
CALL  LNGRIN  ( OFR . YN2 , R , YN2B , D YDX , 6 ) 


LL  LNGRIN  (OFR  , YNO ,K , YNOb > DYDX , 6 ) 
Y ( 1 ) =  YCO 


2 

3 


Y ( 4 ) =  YH20 

Y ( 5  )  =  Y02 

T(6 )  =  YN2 

Y ( 7 )  =  YNO 

SUM  *  0, 

DO  20  1=1,7 

♦ 

xr mri.iT.o.i  iv,to 

19 

Y ( I ) =  0. 

20 

DO  30  I  =  1,  7 

30 

~YTT7"="YTTY"7~5UPT 

RETURN 

* 

TITO . . 

SUBROUTINE  TERR ( X »RS t VD »NSET . RR ) 

UIWEnSTOn'X'  f 2  )  ♦  K ST 40 1 7 VO  ( 4  OT ♦  NS  E'f  72 17  RR  (  2  ) 

LBB  «  NSET ( 2 ) 

- LBJ  =  N'5  E  T 1 1 ) - 

J  *  1 

VS-*-V"j-2 - - - - - - - 

00  20  I=1»LBJ 

— K-v-yrr+K - 

I F ( RS ( N ) -X  (  J  ) ) 20. 30*40 

70“ C'OFJT'I  NUE - 

WRITt  (6*25 )X(J) 

2y‘TrORWfTlH07E"r878r22H__AbOVfc""R'ANGE""0‘F‘‘fABL"t  7 
30  RR ( J  )  =  VD(N) 

'-Tg-J"=T6B - 

j  =  j  +  l 


H-  ( J“3 > 10.30*30 

40  RR ( J  )  =  (X( J)-RS(N-2) ) / ( RS ( N ) -RS ( N-2 ) )  * < VD < N ) -VO ( N-2 ) >  +  Vu(N-c) 

lSJ'V-Ubb  — 

J  =  J  +  1 

TFTj-"37Io’»!To*"5S" 

50  RETURN 

- ENT5 - — 


FUNCTION  XDEL(RVAP) 

- TT1-RV-4p--;T1-r0-r5V5 - - - . - . . . — 

5  XDEL  =  0.005 

- GU  TO  1U0U - 

10  IF(RVAP-. 3)20*15. 15 

- - - - 

GO  TO  1000 

- 2v7-TFTRVAF-F."2'nU'.75"."2'5’ - 

25  XDEL  =  0.02 

- GO  TO  lOoO - 

30  IF(RVAP-0. 1)40. 35*35 

- 3T'x~OFl“'="'0".'o5' - 

GO  TO  1000 

- 4'J-TF1'RV-AF'-0'."051T0"*'4'5';'5'5 - 

45  XDtL  =  0.100 

- GCTTlTTuTC - 

50  IF(RVAP-0.02)60.55.55 

- 5-5"xOFl”=“0".70"0 - 

GO  TO  1000 

- 5v-xD'Et"^”075Tro - 

1000  RETURN 


r 


FUNCTION  PVAP  (  T.  I  ) 

■VAF0Pr‘PRr55URT"aF“NY04TT=TTVMWHfr=2 )  - 

DIMENSION  A ( 2 ) *ti ( 2 ) .C(2) 

■WTTS^D .47  *  14.  3207  >  *  I  b™JT6 'iO . .  7 36 3 . 2 1 )  .  I  C=  1&  1  •’» -&T.  I2T3T 

PVAP  =  EXPH  A  (  I  )-B(l  )/(  r+C(  1  )  )  I 

'RETURN -  -  -  ' 

tND 


IF(UPA)  _9S.  92.  92 

TER'MTN'ATl  ON"  POTnT"“AnO  F6F  WRI  TE  FOR  PLOT  TAPt 


CALL  PLOT  (8.3.  0*  -3) 


FUNCTION  TdL(P.I) 

C  bOI  Ll  N(j"T"£MP'  OF  N2041  f  =T)Y"  'MMh  (  1=2) 

DIMENSION  A( 2  )  »d(2 ) »C(2 ) 

DATA  (A  =  20. "47.T4-.-3"iir.  (d- 12630.  ,  /  363 . 2  >  .  ( t=  1 8  i  .  , -6  3  .  I <<1  ) 
TBL  =  B(  I  ) /(  A(  I  1-LUoHP)  >-C<  I  ) 

RFTORN  . 

END 


FUNCTION  RHO  ( T ,1) 

'IT  I  QUID  WNS  TT  Y”0F""N704TT=  l )  V'Kmh  H='2T" 
DIMENSION  A(2).d(2).C<2) 

DAI  A  ( A=0. U6 16 .0T030623 ) .  (b  =0.1082.0.40289). 
1  (C  =  -0. 5)4321. 0.0) 

‘  RHD '  ="  ‘  A  rrr+H  TTlTf  T71V  E  4 "  VCT I  T#T/Y.  fc~S  *T 
RETURN 

end 


FUNCTION  VISCVIT.I ) 

'VAin5R~VT5CUinT7‘"FUR-"N70'5'TT»Trt  'WWhTT^T) - 

DIMENSION  A I  2 ) »B ( 2 ) 

UA I  A  (A  =  U.4IV98.0. 41998). (d  =  U.OU93bi.O. 00938  IT 
ViSCV  *  -(  A  (  I  )+b(n*T)/l.t7 

'RETURN - 

END 


FUNCTION  CVAP  ( T • 1 ) 

0  VArTOR"5PEXTTTC'UF'N'ZOTf TT=TTV  'MffHlT=27 
OIMtNSION  A( 2 ) *01 2 ) »C ( 2 ) 

- go  ro  i i * i - 

lo  TR=T-1800» 

CvAT5  =  j2-74+T”K*T2V3”u5  t-V-TTR"*5  »  33C.-9T1 
RETURN 


2Tr_rvxpv.-y56+r.'BWt-"4*T 

RETURN 


FUNCTION  FK.A  (  T,  I  ) 

C  VAP‘0ft'"TRET?M4L’~C"ONl5uCTTvTTy"o’F"‘N2O4TrsIT*""MMHTr*’2T 

DIMENSION  A  (  2  )  •  B  (  2  ) 

DATA  I A=-l. 67. 0.12377*  TB='D'."006b2 ,'6. 002,236  > - 

F<A  =  ( A ( I )+6( I )*T ) /1.E7 

RETURN 

END 


FUNCTION  CPL  ( T.I ) 

c  •  ur  ouii3”SPTcrFrc‘'h  t  at*uf"'n2'04  rr=rrv'‘ mmm  <  t=  2") . . . . . . 

GO  TO  ( 10*20) .1 

- 1J  CPL  =0.354  +b. 3E-6  » ( T-430. ) »»? - 

RETURN 

2TJ-CPU"=o;-58_9r25''+ '2Vao'7oft-4"*f  ' 

RETURN 

FND  . ”r  r 


FUNCTION  CD ( RE ) 

r - uRAb"COTr'OT'TJRoP - 

IF  ( RE-80* ) 10* 10*20 

- 10  CD  =27.  *Rt  **(-.57) 

RETURN 

— 2Tj”TFm-i".'r4rjDT3u;T0'' 

3D  CD  =0.271*RE**.217 

- RETURN . . 

4o  CD=2 • 


FUNCTION  ALM  (  T.  I  ) 

C  HEAT  OF  VAPOR  1 ZATTON  UF"N704fT=n'."HM«Tr*2T 
00  TO  ( 10*20 ) *  I 

- TvJ  TR  =  7(T0.~  -1 - 

ALM  =  272. 

IF ( TR  T16 *16* 15  -  - 

15  ALM  =  ALM  +  17‘+.*6uATF  l  Tk/^50.  ) 

lo  RETURN  - - 

20  ALM  =  i jO. 747  - f * ( . J3 V i 505-1 • 2 14E-4*T ) 

Re.  I  URN - - 

END 

-  -  - 7_ 


SUSROUTINE  wTMINT  (  Of-  R  » P  •  wM  .DWMDOFR  ) 

C  MOLECULAR  wT  oi-  CoMdUSTTGN  UAYES  ASU^'EuftCTIUfi'  o?  o/F  "kAT  i u 

DIMENSION  uHAb  <7)*WMTAb(7) 

- D'/n"A  ■  roFT'Ao"= — onr* — rru; — rry; — T7tr, — 7751  jtg  » — 5ToT* — 

1  (wMTAb  =  14. 00*14. 67*19.75*22. 93*23.Bl* 23. 00*23. 00) 

CALL  LNuRIN  (  Ork  *  wM»uFTAb  »  WmTAS  *  Di/MD0F1?*T ) 

RETURN 

- END -  - 


FUNCTION  A ( X ) 

C - d?US3“EErn-0N-AL”A-R-E^-UF-CTfAHBER‘V5"r 

A  =  107.  -j.55*X 

- RETURN - 

•  END 


SUbRGUTlNd  MASStT ( tMO »KHO *K S * EMDOT *N *DkOPN ) 

DTWtH'Sl_oN'-EMoT2U'TVK_sT"20T*L)fioPfrr2“51  7 

CRri0  =  4. 16b7'y0204*RHO 

- rRMTiEHTOT/FLOATF  l  Ml - 

— FRArcT'™F]?/rc/T;"(J - 

NP4  *  N+4  _ _  _ 

EMO ( I  )  =  CRnU*RS( i )**3 

*!  DROPNIU  =  FRAClVtMOm - 

DO  b  1=6. NP4 
EmOTT  )  ="£  Rh  0  *  R  S  IT  )*  *3 
5  DROPN ( I  )  =  FRAC/EMO ( 1 ) 

- RETURN - 7 - 

END 
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bUbKOUT IN t  VUSuhs  ( ur k . P » T « V 1 Sb »  1  Kb » CPu ) 

VISluSITY.  Tr-itki'-iHL  uuNDUC  T  i  V  i  T  Y  •  AND  bPtCJFiC  hlaT  or 
o  o UklbU  5 T  1  ON  OMbc.0  Mb  a  PUNCllON  Oh  0/F  HAliu.  PktbbUkr »  A‘,„ 

c  l  uMPlkA  fokc. 


DJfitNbloN  TTAni  7)  .  v^o(  V  )  .  VCu2  (_7  )  »  Vn*  (  l  )  »  Vi  i2U  (  /  )  ♦  Vue  l  /  )  »VN^  (  /  )  » 
1  VNU  (  /  )  «  5U">  17)*  Wivi  (  / )  *  V  (  7  )  •  Y  (  /  )  .  Vk  (  I  )  *  wivik  (  (  )  ♦  Art  i  (  /  )  * 

£  YPHI ( 7) *oP( 7) 


DATA' 

TTTA6”'= 

50U.U«100u.U»2uGU.O»3000.0* 

4000.0.3000*0. 

7 OoO . 0 ) . 

1 

( vco  = 

0.971.  1.633.  2.51. 

3.34. 

3.90. 

4.31. 

3.30  )  » 

> 

(VC02  = 

0.791.  1.423.  2.32. 

3.2  1. 

3.73. 

4.19. 

4.70)  ♦ 

3 

( VH2  = 

0.464.  0.730.  1.16. 

1.64. 

1.95. 

2.22  » 

C  .  O  3  )  ♦ 

4 

(VN2U  = 

0.312.  l.iuO.  2*22. 

3.34. 

4.3b. 

3./3. 

6 .0  0  )  . 

3 

( VU2  = 

1.100.  1.7b7*  2.63* 

3  o  73  . 

4.43. 

3.13. 

6.30 ) » 

5 

TVNT  = 

U.9ou«  1.3^3*  2.32. 

3.  AT. 

4.71  . 

4  .To . 

5.62). 

l 

(VNO  = 

0 • 03  o  .  1 .  i 1 2  *  2.13* 

3.0^ . 

4.10. 

..30. 

3.30)  . 

tt 

( WM  “ 

26*.  4M.»  2.U16. 

16*. 

32.  . 

2  b  .  * 

30.  ) 

CALL 

LNGRIN 

( T.V( 1  )  .TTAB.VCO.DYDX. 

7) 

TALC 

lNGr 1  ft 

TTTvTD  .TT7TB.vTGT.DYDx 

»T) 

CALL 

LNGRIN 

( T ♦ V ( 3  )  .TTAB.VH2.DYDX. 

7) 

CAUL  LNGRTN 

rr.  VT4T *TT  A6  * VH2U »  UTDX *T7 

CALL 

LNGRIN 

( T .V ( 3  )  .TTAb.V02.DYDX. 

7  ) 

- crnnr 

■■mMnBU-.n'nwwin.f 

l  J 

CALL 

LNGRIN 

( T ♦ V ( 7  )  .TTAB.VNO.DYDX. 

7) 

TAL'C"  TUTCFtTC>FF;f;Y7 

SUMVIS  =  D. 


If-'  (  Y(  1  )  .LT.O.  001  )  20*  9  » 

~  buRTTl  =  u.  T 

DU  10  J=  l » 7 

Vfi-rjr-r  vTff/"vTjT 
WMK  (J)  =  WM(J)  /  ^(1) 

phtu r”"*rr.  +  ssrtttvto  >  *"3ci?TrnTOirrjrrrr**5' 7Ci;eTo~~*~sut<  w  < 

1  1.  +  1.  /  WMR{ J) ) ) 

YPHI(J)  =  Y ( J )  *  Phi  (  J ) /Y ( I  )  , 

lu  SUM ( I )  =  SUM  ( I )  +  YPHI ( J ) 

TumvTs  ”=”  "s'umVIs  +“  vTrT"/”s"uMTIT 

2u  CUNTINUt 

VlbG  =  SUMvTs  /‘iitb 

C P ( 1 ) —  9.46  -  3290. /T  +  1.07t6/T/T 

CP ( 2 ) =  16.2  -  0530.  /T  +  1.41L6/T/T 
CP ( 3 ) =  3.76  +  3.7  6  E-4*T  +  20./SQRTF(T) 

”  cp ( 4)=  ir.'b6'-“v^:7iwrm”75oor7"T'' 

CP  (  5  )  =  11.313  -  _  _  _  173./  SbRTF  (  T  )  +_  1530./J 

■CP"(6)"=”6."9362'  +~ '3.~3“7 3b-"4"'*"f 

C  P  (  7  )  =  7.7203  +  2.463D-4  »  T _ _ 

SIGA  =  0. 

SIGB=  u. 

D0"3u  I  -  j~»  7 

SIGA  =  SIbA  +  _  Y(  I  )  _*  C_P  (  I  )  _ 

TJ'  *  sTgb  "=”sTCtT””+"  y  fi  r  *  wMTTI 

CPG  =  SIGA  /  SIG6 

- “ITu — =  VibG  *  I  CPG  + — 2.4b'T  SIGd) - 

ktTURN 

- E-NU - 
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SUBROUTINE  LNGR1N  (  XARG*  YARG*XTAtJ*  YTAB*DYDX*N  ) 

”  ”‘CT0X0lfATir‘LAN(3ft^N(SlWlWRF0UArf0N"0ft"ExfW0LXfrON'(fOR  F  ( x j 
AND  OF ( X ) DX 


Irtdl30l  *  Y  I  Ad  ( 

IF(XARG.LT.XTAd(l)  )  10.  20 

17  "  WR l‘T E" ’  T6 *"f  015 1 " ~X  A RG  .  . 

10U  FORMAT  (10H0  XARG  =  F12.6*  41H  lb  OUTSIDE  THE  RAN  OF  TAdULATt 

-----IIJ-VALWS;7By;7lH  VAru¥S"0r-rryrAND"r'TxT_TfETuTFNrD'W‘rALLlNG  TPROo 
2RAM  ARE  EXTRAPOLATED.) 


GO  TO  50 

'27  TF  "  f  5CA  ft'G  .7  TV  x'  TA  oTNl  T  ”  "i"  0  V 4u 

3d  WRITE  (6*100)  XARo 


-3T"”  r-*—Fr-"r- 

GO  TO  50 


f:cr<i;n 


41  IF  (XARG. LT.XTAd(N) .AND. XARG. GT.XTAd(N-l ) )  31»42 


— M--£-7i”-"r - 

DO  43  I  =  3*  M 
'TF'"naA'Drrr.‘GT;7AT<GT‘757"'zr3' 
CONTINUE 


n:ici: 


DIF2  =  ABSF  (XARG  -  XTABU  +  1M 

•TF"rT5TFr".Tr;DTF7T""4T;"5F - 

K  =  I  -  2 

'TTO_7Cr50 - 

K  =  I  —  1 


B  *  XARG  -  XTAB  if.  1 ) 

■T'V-^AirG"-_7T/T57R+2T - 

R  *  XTAB(K)  -  XTAtt(K+l) 

“T’vm'BTKy'-'xTwrK+TF - 

T  =  XTAB  (  K  +  l )  -  XTAB  (  K+2  ) 


R  F  (  X  ) 

YARG  =  YTAB(K)  *  d  *  C  /R/S 


*  A  *  B  /S/T 
DIFFERENTIATION  FOR  b  (X) 


-  Y TAB ( K+l )  *  A  *  C  /R/T  +  YTAb(K+2) 


^T)X7'7T7B7M“7'T6VCT7R7s"-'’rfAliT'K+Tr'*""rA"+c7/'R/T'+"YTA"B“('x+2y 

1*  (A+BI/S/T 


irfi-ic 


I  *  m  #■  l  :VB  an*] 


60  WRITE  (6.200)  YARG*  DYDX 

'Ton  ro^A7-727x7?HYAft7-v-rr27Fnox7rHDYDx'»'"r2oriov 

61  CONTINUE 

TTETOITN  ‘ 

END 


SUBROUT  I Nfc  TG1NT  ( Or R , PC t T » uTDOFR ) 

C  ADIABATIC  FUAMt  TEMPERATURE  AS  ‘A  FUNCTION  or  u/r  raTIu  t 
C  PRESSURE  EFFECTS  ASSUMED  NuG*  FOR  FIRST  APPKuX.  (  Iim^LululU  l*Ilk> 

- D  rWElTSTT UTTTCTTS  T  *  "UFTAb  I  yj - 

DAT  A ( OF  TAB=  0  •  0  *  0*94*  1 • U  *  1 • 2  *  1*4*  1 • 3  »  i  •  /  5  »  ^  •  0  »  > •  u  )  * 

I  (TCI  =  4T60  •  »  47E0 •  *50B2  •  »5 51  b  •  *  565 <5»  *  5b50«  » 5  jV 3  •  *  :?*:o  7  •  *  jo c  u  •  ) 
CALL  LNGRIN  ( OFR • T »UF TAb . TC 1 *DTD0FR *  9 ) 

"  "'RETURN - - -  ' 

END 


FUNCTION  RT(RS) 

- C - CALT:UL7rrES"'DT!>rAISrCE”OF"UI3'SU,CrATrON"TC.7fHE"WTUNCTTON"aF'  k  AD  I  oS 

DIMENSION  RRS ( 5 ) • RR 1(5) 


- DAI  A  l  RMS  = - U.C1 .  iUuTWuu.  »  joU.  OUO.  )  » 

i  ( RRT  =  HiOt  13»*  i  7  •  t  16»»  lti«) 

C  CONVERT  R"S1  iCJ  TJ  MICkURS" 

RST  =  RS 

TTST  =R3TT 

CALL  LNGRIN  ( RST » R T * RRS »RR T » DRTDRS * 5 ) 

k 1  =  RT  /  2*S4t+4 

RETURN 

- - - ENC - 

FUNCTION  XK ( Ro ) 

c  calculates' dTSsucTaTIUn  Burning" rate  constant  As' function  wr 

C  RADIUS 

- orHENS10(5-  XXKI5)  ,RRS(5) - : - - 

DATA  (XXK  =  3*2*  b«3»  7»3»  ti»2»  11»)» 

•  IRS  =  0 •  0  » 1 00 • * 2 00 • » 300 • » 300 • ) 

C  CONVERT  RS(X)  TO  MICRONS 

"  'RsT'v'rs  . . 

RST  =  RST  *  2.S4E+4 

Call  lngrin  (rot  *xk»rro»xxk»uxk.dR;>.5  ) 

XK  =  XX  *l.E-3/2. S4/2.54 

RETURN 

END 


function  mP(x) 
ap  =  -3'.5V 
RETURN 

end 
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SUBROUTINE  LOGNRM ( RM»SlG*N*RS) 

c - coMPirrrs"RADrr"or^rr‘crROP’S“'ABouT"WFfrcH'T7'N‘Tfl'‘or _THr"FAS3 '  ts 

C  CENTERED  IN  THE  LOGAR I THM I  CO-NORMAL  DISTRIBUTION  WITH  MASS-MtAN 

C  RADIuS  Rm  ANu  SIANDAKU  DEVIAIIUN  "LN  S1G.  AM*  I  HE  INVtKSt 

C  PROBABILITY  FUNCTION  IS  APPROXIMATED  BY  A  RATIONAL  FRACTION. 

- inwE7rsTON"RrsTTcrr - - - - - 

D£L=1.00/FLOATF(N) 

- -STGUN’SL-OGFTSTGI - - - 

P*-DEL/2 • 

- rcp^.— ^ - 

MED  =  (N+D/2  +  4 

- -jv-Tjpzr - 

DO  4  I  =  1.NP4 

- TTTr^rZritrnT--TI - 

10  P=P  +  0.6#DEL 


uO TO  2 

11  I F ( i -5 )  15*15.20 

15TiFfTT,T*r)EL‘ - 

GO  TO  2 

'2TJ-FiFfUET - 

30  I F ( I *GT .MED) 5*2 

"  Z  1T5UK  rFI  -LUGF  <  )  ) - ’ 

XP*T-<  2.3075 3  +  . 27061 *T)/ < 1 • +  T* ( • 99229+.044ttl*T ) ) 

- rSTG3riETFFT3TGLTT»TPT - 1 

1FC 1-3)4*3*40 

WTrrr-"5T4-,-4Tf5TJ - 

49  P  *  DEL/2.0 

- GQ-TCJ  "g - 

50  J«N  -  I  +  9 

3'R'SU  r«R~M*ES"fGX  . 

4  RS( I )«RM/ESIGX  _ 

"5""R"ETUR'N  .  . 

END  i 


SUBROUTINE  DIFt-U  (UFK.P.T  *VPU2*i  »UFU) 


c  urrFu5r^n?"'DT‘NrjMiviT"V'wrr=2)"TN"'vAPoR''F"rLM  . 

DIMENSION  FACAI7) .FACB (7) • Y I  7 ) • EPSA ( 7 ) . EPSB < 7 ) »ETAB( 14) . 

- 1 — oMbsnwr  pm  7) - 

DATA  (FACA  =  4.B6*  3.o7.  i/.D.  7 . 4 i .  5.b:>*  4 •//*  4.Yj). 

1  rEPYA  =  344  •»  4Y3  •  *  2173  .  •  J » »  J4V. »  VIV. »  3  3  o  »~)  * 

2  (FACb  =  4*77*  3«ou»  17.2*  7.27*  5.7a*  4.o/»  4.o4)» 

3  - rEPr5"'="4zi;T"b"5T;v'Y87r;y"Td5vv'4irav'."'3"tt6vr'43'd"."rv'“ 

4  (ETAB  =  U • 3  *  1.0*  2.0*  3.0*  4.0*  6.0*  b.O* 

5  ro.o,  20.0.  3u.U*  4(5.0*  60. u*  ttO.U.lQO.Gl* 

6  ( OMEGT  =  2. 622  *1. 439*1. 075*0. 949* 0.864.0.612*0. 771. 

- 7 - Tj;T47^;¥47r;7r;v2T?rr:5^670V5y9V0V53"5V0V5T7r"' 

CALL  COMP  (OFR.P.Y) 


— UO"TO“‘U*T*T - 

10  Y ( J )  =  Y ( J )  *  (P-VP02)/P 

— sumo  «  j; - 

IF  II.EQ.l)  20*40 

2V-y«UW-~i;-~VTO77F - 

DO  30  J  =  1 » 7 

— r-v-T_T-rp5-ArJ1 - 

CALL  LNGRIN  ( E .OMEGA , E TAB .OMEGT *DYDX » 14 ) 


FUIJJ  a  r  ALA ( J  J  *T  **1.5  /OMEGA 
30  SUMD  =  SUMD  +Y ( J )  /PD(J) 

40  XNUM  *1.  -  VP02/P 

■ — TTO"5Tr;r=T."3r - 

E  =  T  /EPSB(J) 

- CALL'  LhGKTKI . I E  *  O'M'  L  G  A".  E  7  A'dTSWt  oT*  JYT)  XTITF 

PD ( J )  *  FACb ( J ) *T**1 .3  /OMEGA 

'5o*‘"S0ffD'-a"“5u"Mu""V'yTjT/TroTj') 

60  DFU  =  XNUM  /  SUMD  /P  /l.c6  *  14.7 

— ‘K’ETuEn - 

END 


S .  Sample  Problem  with  Output 

A  sample  problem  was  solved  and  typical  program  output 
is  presented  here.  The  input  parameters  are  shown  in 
Figures  20  and  21*  The  following  pages  show  the  printout 
of  droplet  histories  at  arbitrary  x  increments  along  the 
chamber  length.  Printed  output  in  this  solution  involve: 

(a)  State  of  x  position  (inches). 

(b)  Calculated  gas  velocity  u, 

(c)  Chamber  temperature,  T,  and  pressure  P. 

(d)  Calculated  local  bulk  gas  O/F  ratio. 

(e)  Local  Mach  number  of  gas. 

(f)  Cumulated  fraction  of  0,  F  and  both  vaporized. 

(g)  Local  vaporization  rate. 

(h)  Drop  size  and  mass  for  each  group. 

(i)  Drop  absolute  and  relative  velocity. 

(j)  Droplet  liquid  temperature. 

(k)  Reynolds  number  of  each  drop  based  on  the  speed 
of  sound. 

Figures  23,  24  and  25  are  samples  of  the  printed  output 
at  various  chamber  stations. 

In  addition  to  the  printed  output  the  parameters  of 
most  interest  are  machine  plotted.  These  plots  display 
trends  in  these  parameters  in  a  more  meaningful  format 
than  the  printing.  The  results  for  the  sample  case  follow 
the  sample  printouts. 
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APPENDIX  II. 


COMBUSTION  INSTABILITY  PROGRAM  DESCRIPTION 


The  Dynamic  Science  Corporation  combustion  instability 
program  solves  the  following  set  of  nonlinear  partial  dif¬ 
ferential  equations.  Cylindrical  coordinates  are  used,  in 
which  the  z  direction  corresponds  to  the  axial  direction  in 
the  chamber.  The  following  dimensionless  equations  are  the 
result: 

Continuity : 


(ll-l) 


Momentum  (in  8-direction): 

3v 


•  Uli  ♦  ’’'I9  .  1  »>"  „  L«'v',fh)  •  «  i>><' 

9t  e  Y  98*  ^  (96’) 


6  f(Y) 
2 


f I 1-2) 


Ener^n 


’’($■  -J  'y*  nf  (t*u  '■(&  *  £)• J  tSV1*’ 

♦  L««  J(Y  -T')  ♦  y[(Av')2  ♦  V'e2]j  f(Y) 


CII-3) 


Ideal  Gas: 


9P' 

TF*  *  T  Tf*" 


9T» 

♦  #'  TT 


(H-4) 


The  set  of  equations  (II-l,  2,  3)  represents  a  mathe¬ 
matical  model  used  for  the  analytical  determination  of  the 
minimum  pressure  or  velocity  perturbation  required  to 
develop  into  a  traveling  wave  within  the  combustion  chamber. 
The  one  dimensional  model  employs  an  annular  section  of  the 
combustion  chamber  of  thickness  Ar  and  of  length  As* 
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Previous  solutions  to  analyze  combustion  instability 
employ  a  first  order  explicit  finite  difference  scheme  in 
the  time  direction.  The  spatial  or  theta  derivatives  are 
approximated  with  a  Stirling  first  order  central  difference 
relation.  Initial  investigation  demonstrated  that  indications 
of  combustion  instability  were  not  entirely  valid.  That  is 
error  propagation  during  the  computer  solution  was  giving 
false  indications  of  combustion  instability. 

Since  the  solution  of  the  system  of  second  order  non¬ 
linear  equations  essentially  represents  an  initial  value 
problem  for  a  system  of  differential  equations,  it  was  decided 
to  use  more  stable  but  numerically  more  cumbersome,  higher 
order  difference  methods.  The  Dynamic  Science  computer  program 
uses  a  third  order  Adams  Bashforth  predictor  given  by: 


yk*l  "  yk  ♦  yj  (Ssfk-s9fk-l°7fk-2-9fk-3)  (II-S) 

and  a  third  order  Adams  Moulton  corrector  formula  of  the 
form 

Xk*l  *  Xk  ♦  TT  (9fk*l+l9fk-5fk-l*fk-2)  (II“6) 

to  advance  the  solution  at  each  time  step.  It  should  be 
noted  that  previous  methods  of  solution  to  the  nonlinear 
set  of  equations  involve  only  first  order  methods  to  advance 
the  solution  to  the  next  time  step.  The  present  difference 
equations  require  a  history  of  four  previous  points  in  order 
to  advance  the  solutions  to  the  next  time  step.  The  first 
order  method  sustains  an  error  term  of  order  (At)2  while  the 
predictor  corrector  formulae  (II-S  and  6)  sustain  an  error  of 
order  (At)5  and  hence  allow  larger  step  sizes  during  integration 
while  maintaining  numerically  stable  solutions  and  minimizing 
error  propagation. 

The  Dynamic  Science  Corporation  program  corrects  only 
once  and  rather  than  continuing  to  iterate  to  converge  to  the 
solution  a  test  is  made  during  the  integration  cycle  to  check 
for  four  significant  figure  agreement  between  the  predicted 
values  and  the  corrected  values.  If  sufficient  agreement 
exists  then  the  solution  can  continue.  If  sufficient  agree¬ 
ment  does  not  exist  then  the  time  step  is  reduced  by  a  factor 
of  2  and  the  solution  is  continued  from  the  last  point  which 
satisfied  the  four  significant  figure  agreement.  If  sufficient 
agreement  is  maintair.ed  for  eight  or  more  points  provisions  are 
made  for  doubling  the  independent  step.  At  the  present  time 
the  step  control  la  based  on  the  absolute  difference  between 
predicted  and  corrected  ordinates  of  the  integrated  functions 
p,  v,  and  T, 
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In  tha  spatial  or  thota  direction  tha  accuracy  of  tha 
that*  darivativa  has  baan  iaprovad  bv  aaploying  highar  ordar 
cantral  diffaranca  foraulaa.  Tha  first  ordar  prograa  usas  tha 
following  aquations: 


(jT)  Ta*l tn"Ta- 1 ,  n 

TaeTa^n  " 


(II-7) 


Whila  tha  Dynaaic  Scianca  Corporation  coaputar  solution 
aaploys  tha  following  fourth  ordar  cantral  diffaranca  rala- 
tion. 


TOT 


Ta42,n| 


(11*8) 


whara  T 

a(n 


T 


(a&e.nht)  * 


Tha  darivativas  in  tha  axial  diraction  ara  determined 
with  tha  assuaption  that  tha  total  nass,  aoaantua,  and  anargy 
in  tha  annulus  ara  constant.  Furtheraore,  it  is  assuaad 
that  thasa  derivativas  ara  indapandant  of  r  and  6.  Thasa 
assuaptions  rasult  in  tha  following  aquations,  which  parnit 
avaluation  of  the  darivativas  taken  with  raspact  to  z: 

Continuity : 


US 


N 


Energy : 


Equation  (11-10  -  12)  lead  to  the  following  system  of 
algebraic  equations  to  be  solved  at  each  step  in  time. 


*lxl 

♦  a2x2 

•  ci 

•3*1 

♦  »4*3 

2 

•  c2  ♦  ajxj 

(U-13) 

*6  *1 

♦  1 L  *2 
y 

♦  1*  X, 

7”  *3  -  c3 

87*2  ♦  *8X3  ♦  2wx4  ■  0, 

where  (xl#  x2.  x3,  x4)  represents  (  |^.,  |£.,  |I,  ) respective iy 

and  •} • *  *  * • *8  «re  obtained  from  (11-9,10,11,12).  It  might  be  noted 
that  the  z  derivatives  are  permitted  to  chang  :<  with  time  but 
remain  constant  around  the  annulus  at  each  time  step.  (i.e. 
only  gross  property  changes  are  permitted  in  the  solution  of 
this  mathematical  model.) 

The  integration  cycle  consists  of  six  subroutines  defined 
as  follows: 

(1)  Subroutine  ASET  determines  the  coefficients  of 

system  of  equations  (J-13  to  obtain  the  axial  derivative-  based  on 
the  assumptions  that  the  total  mass,  energy,  and  momentum 
remain  constant  in  the  annulus  and  the  axial  derivatives  are 
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independent  of  r  and  6. 


(2)  Subroutine  ZDIR  solves  a  fourth  order  nonlinear 
s/stem  of  equations  for  the  axial  derivatives  using  a 
second  order  Taylor  expansion. 

(3)  Subroutine  TDIR  computes  the  t  derivatives  at  each 
of  the  theta  nodal  positions  using  equations(II- 1,  2,  and  3). 

(4)  Subroutine  NADM  uses  a  third  order  multistep 
procedure,  equations  (H*5  and  6}.  to  integrate  the  variables 
p,  T,  and  v,  and  compute  the  absolute  difference  between  the 
predicted  and  corrected  values  for  automatic  step  control  at 
each  of  the  theta  nodes  of  the  annulus. 

(5)  Subroutine  TUP RED  determines  theta  derivatives  of 
the  integrated  variables  for  each  of  the  theta  positions 
around  the  annulus  using  equation  (II-8). 

(6)  Finally  subroutine  SHIFT  moves  the  computed  points 
at  the  current  time  step  to  the  position  of  the  previous 
time  step  in  order  to  begin  the  integration  cycle  at  step  (1) 
for  the  succeeding  time  step. 

The  remaining  subroutines  can  be  described  as  follows: 

(1)  REED  -  reads  initial  data. 

(2)  ORG  -  initializes  print  and  integration  counters. 

(3)  RSET  -  Initializes  pressure  velocity,  temperature 
and  density  around  the  annulus. 

(4)  AVGE  -  Forms  (PMAX-PMIN) /PAVGE  for  plotting. 

(5)  DRAW  -  Subroutine  to  fora  graphical  output. 

Figure  38  represents  a  graphical  presentation  of  the 
logic  involved  in  the  main  program.  The  routines  involved 
in  the  integration  loop  are  described  above.  The  test  on 
MPTN  indicated  in  the  previous  diagram  is  merely  a  print 
indicator  in  order  that  the  print  will  occur  at  equal  intervals 
even  though  the  independent  step  may  be  changing  in  magnitude. 

Figure  39  represents  the  input  card  format  for  the 
Dynamic  Science  Corporation  combustion  instability  program. 

The  graphical  display  of  the  indicated  input  data  is  shown 
in  Figures  40-49.  Finally,  Section  II  contains  a  complete 
listing  of  the  Dynamic  Science  instability  program. 
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* 


FORTRAN  CODING  FORM 


card  form,  IBM  electro  888157,  is  available  for  punching  so  wee  statements  bom  this  form. 


FIGURE  40 
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PRESSURE  WAVE  PERTURB 
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FIGURE  46 


PRESSURE  WOVE  PERTURB 
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PRESSURE  WAVE  PERTURB 
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MAIN  PROGRAM  FOR  COMBUSTION  INSTABILITY  MODEL  OF  DSC  _  _ 

BTmTnSTon  f c7¥olT NMtTdrr" bd ( io*2l.i6 ).  art”  zooo»2 ) *  v  ( 1  o  *  z  a  ) 

1  ♦  P  ( 10*21 )  *RHO  (10*21)  *  T  <10.  21)  _  _ _ 


10*21 ) *DTTH< 10.21) .ORhOTHi 10.21 ) .DVT (10.21) 
1  iDVTH ( 10.21)*  DRHOT (10*21) 


DIMENSION  W(  10  *21)  *  WZdO.21) 

^COMMON  / E/  FC.NM  /B/  BD  /C/  ART  __  _ 

E0UmrENCE"(trDT294iT,  W)T  <6 573 151” *  WZ) 

EQUIVALENCE  (BP(1051).DTTH)*(BD(1891)  *DRH0TH)>  (BPtZSZl )  »1>VTh 


EQUIVALENCE  <FC(1).  TI).  <FC<31)*  TSTOP),  (NM<4).  MGAM )  » 

_ 1  ( F C  (  3 )  >  H  )  *  ( NM  (  5  )  »  MPTN1*  <NM(3[*  _MALP  )j_  (NM(13).  lj  »  _ 

2  I F  C  (  3  8  )  •  Z  IP  )T7fT(  32 )  7'x L IVTfC  <  34  )’•  7  J  )  »  "(~F  cl  35  f  .  DfcLV  7 . 

3  ( FC  (  36  )  »  GAM)  *(NM(7]*NSW)  __  _ _ 

equivalence  TboTssTIV  "nTTB'5”8  4T)7”d T f )  *""(B"J(2lbn  »  V)  . 

1  (30(  2311  )  »  DVT)*  ( BD  < 1471 i »  RHO)*  <BD(1681)*  DRHOT > »  (bO<l).  P) 


FORM  INITIALIZES  COUNTERS  FOR  INTEGRATION  ROUTINE 

Li^’jOJj:^.fi2R..E^QmW&.&QUIJNL5./JA_IIIL£§. _ 

EXECUTE  DATA  READ  FOR  BOTH  INTEGRATION  AND  STABILITY  INITlAt-IZAT 
I  *  1 


CALL  REED 
CALL  RSET 


SUBROUTINE  ORG  INITIALIZES  PRINT  COUNTtRS  AND  SETS  UP  THE  NcCESS 
ARY  INTEGRATION  TERMS 


CALL  ORG 

SET  UP  COEFFICIENTS  TO  SOLVE  FOR  Z  DERIVATIVES 

TCf — rarc-ASET - - - 

SUBROUTINE  ZDIR  SOLVES  FOR  Z  OtRIVATIVtO 


CALL  ZDIR 

NOW  SOLVE  FOR  T  DERIVATIVES 


CALL  TDIR 

NADM  _PERFORMS  THE  ACTUAL  NUMERICAL  I NTtGR_AT ION  „ 

---  l  l“na"dmTdr  ho  tT”r"h"6V~T7  “  ‘  ’ 

CALL  NADM  (DV 


CALL  NADM  ( DTT  »  T.  3) 

_ „CAkW__IJiPREP_ _ 

NEXT  TEST  FOR  PRINT  POINT 

_ _ _ _ 

BRANCH  TO  50  IMPLIES  PRINT  POINT  OBTAINED 


Ml 


ra<ia 


51  WRITE  <6*21 )  (ART  <  1  .  J)  .  J  =  l*2) »H. ( P( 1. J ) »RHO (1»J)»T(1.J)»V(1»J). 

_ A.XillP_lj.WZ_L^aLi__J=l^JLl _ _ _ 

21  FORMAT  (1HO*3E20»8/(6E18*8)  ) 

_ ilSXI_E_tiu4^J _ LPJlIllj-OJ.lP.IJHU^lJL.P.eHOJ.LXtJ-LjpHtLOJJjjAjOJj*. _ 

1  DVT ( 1 »  J).  DVTH ( 1  *  J ) »  0*1.21) 


52  I  =  1+1 

_ IiSX_f_QIL...IXML-SiOP.. _ 

IF  (TI  -  TSTOP)  60*  55.  55 

_ SXQilE^RL£LI_J^QlDLX-f.Oil _ £LUJ_lJJiG_ _ 

55  ZIP  =  ART  11*2) 


CALL  PLOTTING  ROUTINES 
GO  fb"~5" 

SUBROUT J_NE  SHIF_T_  ^PDAI^S  TERM.&  INVOLVED  WITH  If 

6b  "call  sh  Tft  ,,n 

_ GO  TO  10 _ _ 

END 


SUB  ROUT  InT'REED"  . 

DIMENSION  FC { 60 ) *  NM(20)»  BD ( 10 » 8 1 » 15 )  ,  WT(20) 

COMMON'  /E/F  C  .  NM " 7 B  /  ri  D . 

EQUIVALENCE  IFC11)>T),  ( FC ( 6 ) »HMAX ) t ( FC ( 6 ) . EMI N )  > _ 

1  ( F C ( 9 ) » tMAX )  ,  ( FC ( 1 1  )  »  WT { 1 ) ) *  (PC<31)»  TSTOP ) TTFcTTZl »Xl) ♦ 

2  (  FC  (  33  )  *  RED  )  ,  (FC(34),XJ),  (  FC  (  36  )  »OtLV)  .  (  FC  (  36  )  .OAM  )  , 

”  ”3 "  Tfc  <  3"  7"f  i  a  p  V r  7  f  c '( '3  9)  \  s’  c')T  ’  (  f"c74T )’  V  ‘  v’ 1 1 1  f"  c'(Ta  )  To  t  h  ] 

EQUIVALENCE  (NM(7).NSW).  (NM(9).MP).  ( NM < 1 1 ) ,N0 ) . < NM ( 14  )  ,  imO  ) _ 

. "d  6"  I  o'  ’  T "  =" "  Y ."  2 6 . . . . 

10  WT ( I )  =  1.0 _ 

READ  (  3.20)  T.HMAX.  EMIN,  EMAX  ,  XL  » RED  »  X J ♦ DELV » OAM,  AP  .SC.VA.OTri 

_ 1..TST0P  _  _ 

TF"(  EOF  ,  6  oT  9  0  *~1 8~~  '  . .  . . . 

20  FORMAT { 6l  12 • ti  ) 

T8"“KrA5T_V;7IT“Nsw%“MP7"N6V“N5'  . 

21  FORMAT  (4112) _ _ 

XND  =  ND  -  1 

_ _DTH__=__6_.2e_3 165  3071  /XND 

wr  iTeTsTaoT  t“  hma x  Vi  mTnVe  m  ax'^x’l’.'r  e"l57x  j75" el" v  t  oamVapTs’c”,  vz7u~f  n’»""" 

_1  TSTOP 

40  "" "for’mat TIriT»"2 IW t7Hy.A7rE"MTN7xL‘rRTD7rj/Th’“r7“£'i'r.'7 / 

1  1H  .26HDELV. GAM .AP.SC.VZ.DTH, TSOP/ 1H  ,7tl3.7) 

rETUrn 

9  0_  END  F I  LE__  1 7 

sTop"  . 

END 


SUBROUTINE  RSET 

DIMENSION  P  (  1 U .  8 1 )  •  DP  T H ( 1 0  *  8 1 )  »  KHOJ.l  Q.,  8.1  J_._  __D  I  T.H.t.lp.*  tt.lj.t _ 

1  DRHOT  (  lOtTliV  D2T"frt(Td*8il  *  "60(10*81*15')*  T(l"6*dl)» 

2  V  (  10 »  8 1  )  *  DVTH1 1 0 » 8 1 )  *  D2 V TH ( 10 * 8 1 )  »  FC(8U)*  NM(20) _ 

3  .  DRHOTH( 1 0  » 6 1 ) ♦  A(8) 

COMMON  /£/  FC.  NM  /8/  BD 

equTv"ale’nc£"”<  tfoTiV*’ p‘i"*" ' ”Tbo ( ¥11" ) »  dpthT*  i  bdT  16  2 1 )  *  t  > *" 

1  ( 6D( 3241 ) »  DTTH),  (80(4051)*  D2TTH ) »  (BD(4861)»  RHO) , 

2  Tbd'I 6'48T)T'DKH0fH  rV”iycT54  )V’b“fN‘)  VTFCT3"6l’»  (jAMTV . 

3  ( BD( 7291 ) *  V)»  ( bO ( 89 1 1 ) *  DVTH ) »  (aD(9721)*  02 VTH ) «  (FC(37)*  AP )  ' 
4. ( FC(59 ) »SRD)  ♦  ( FC ( 60 ) * SCd ) .  ( FC ( 61 ) »SCR ) .  ( FC ( 35  )  * DtL V ) * 

5  (  FC  ( 40  )  *  DEL2  V  )  *  (FC(64].bBl»  (FC(65j*  §CJ  »  J  FC  (  66]  »C  2  j  . 

6  ( FC ( 33 ) *  RED)*  ( F C ( 3 9 ) »  SC)*  ( F C  T 3 2 ) »  XL)*  (FC(34)»  X*>)» 

_ l..Jf.C14A)j„fGAMJ.._.,jFC_(.61lA§J.PltJ-FC(6.3J_t5J.e^lj„J.FS:.L^i.tijLejL _ _ 

B  *  ( FC ( 43 ) *  A)*  ( F  C ( 4 1 )  *  Vi) 

_ EQUIVALENCE  (NM(14)*  ND).  ( NM ( 15 ) »N1 ) » ( NM ( 20 ) >NOD ) »  (FC(70)*U2) 

I*  ( FC  (  7 1 ) »  DSQ) ,  (NM( 16) »N2) . £NM( 17) »N2) »  (NM(l6)»N3) 

2  •  <NM(19)*  N4 )  *  (  NM  (  12  )  *  N5  ) _ _ 

c  7 

C  SET  UP  INITIAL  ARRAY  _ 

. 

GAM1  =  GAM  +  1.0 

FGAM  *  SQRTF ( ( 2 • 0/GAM1 ) ** ( GAM1/ ( GAM-1.0 ) ) ) 

BB  =  XL*FGAM  _ _ _ _  _ _ 

bZ  =  1.3333333333*BC 

-----  —  - 

_ SIP  =  SIP2/2.0 _ 

ZIP  -  GAM  +  SI P*D£L2V 

_ SRD„-„SQRTF.L  RED;. . . . . . . 

SCB  =  .6  ~  *SC**~.  333 33 3333  33 

_ _ _ S_CR_  =  S.CB*  SQRTF  (  DEL  V  ]*SR_D  +  £.  0  _  _  . . . 

A  (  2  )  =  6.28316b  3  0^7  1  * VZ 

_ Alt)  s  f.SAM _ 

D2  =  12.0*DTH 

... _ _ 

NOD  «  10*ND 

_ _ 

NZ  =  ND  +  3 

_ N2  =  NZ  +  1 _ , _ 

N3  =  N2  +  1 

_ _ 

N5  =  N4  -  1 

_ £_Qtil„= _ ljLOy.QAM _ 

C0N2  =  1.0  -  CONI 

_ DO  20  I  =  1 »ND _ 

XI  =  I  -1 

_ .Z.I.G _ =  YI»DTH _ . _ ^ _ 

P ( 1  *  1  )  »  AP*SINFLziG)"  +  1.0 

_ _DP.IHLl.ti  L__=„AP*Cp.Sf  LZ  I_s_ l _ _ 

T(i.f)  =  P  ( I  *  I  )#*C0_N2 

_ RHO ( 1 » I ) =  P ( 1 » I )»»C0N1 _ 

V( 1  *  I  )  *  U.O 


m 


DIMENSION  FC ( 80 ) »  NM ( 2 0 J *  W(20) 

(T6MMON-'7r7‘"FC""NM  . 

EQUIVALENCE  (FUII.TIi  (FC(2).RKT).  (FCt3)»H)»  (FC(4).HO)» 


(  F  C  {  P  )  ,  HMI N  )  »  (FC(6!  ,HMAX)  »  (  F  C  (  7  )  »  »L02  )  *  (  FC  (  b  )  » trt  t  N  )  , 

2  ( <*C  ( 9  )  »£MAX  )  *  (  FC  (  1 1 )  #  W  )  ♦  (  F  C  (  10  )  *  E 1 )  * 

TTN^riTrrMiv“ii\M7rrviTDKTr'(N'MT3T*PrALPrrTNMr4r.Mu'AMT*'" . . . 

4  ( NM ( 5 ) »MPTN ) *  ( NM ( 6 ) ,MPTS ) *  (NM(7)»NSW)»  ( NM ( 8 ) »NCOU ) • 


5  (NM(9 ) »MP ) »  (  NM  (10)  »NV  )  ♦  (NMU1J*  NO) 

**##**********#***#*jHHH(*t»fr#**#*#*#***#t**##*##**###****#*#***#< 


DESCRIPTION  OF  THE  LISTED  VARIABLES 
T  r_THIS  CELL  CONTAINS  CURRENT  INTEGRATION  TIMc 
"T<T  ”"ST~ ART "~TTme“0R  ' " PRE'vT 6u's"'rUnge' - <U f TA-  T TmE 
H  _T  .CURRENTLY  USED  STlP  SIZE  IN  COMPUTING 

HO  -  STORE’D  ST_EP"”sTz"£  . 

HZD2  -  HALF  OF  STORtD  STEP  SIZE 
HMI N  -  MINIMUM  STEP  SIZE 
HMAX  -  MAXIMUM  ALLOWABLE  STcP 
EMI  N-WMXr-M7N"77"MAx"’A"LL0wAerE"’ER"R0R 
W  -  ARRAY  OF.WEIGhTS  TO  wtIGH  C.RROR  CONSI  DtRAT  I  ON 
’"”rM‘"-"7o‘7r"76oL5’7o"rN'fs’"FK"oy7"R7<‘‘s"f  art  . 

INDR  -  INDICATOR  FOR  C.RROR  OUTSIDE  OR  WITHIN  MI  N  MAX  TulEkANCl 


-  COUNTER  FOR  RU.NGt  KUTTA  I  NT  ERMED I  ATc  POINTS 

-  PHASE  INDICATOR  - 1 , PRED I CTORO ♦ R . <  1»  CuRReCTuR 


N  -  PRINT  COUNTER,  CURRcNT 

S  -  TOTAL  NO  OF  POINTS  IN  PRINT  INTERVAL 

"”-p-R7Nf-TN(51i:A7oR‘T77Ao7M"R0u'frNE' 

U  -  TOTAL  NO  OF  COMPUTED  POINTS  DURING  I  NT  UOKAT I  ON  CYCLc. 


INITIAIZED  BY  AN  ADDITIONAL  ROUTINE 

it******************  *  *****  **«•  ************************************* 


HMI N  =  HMAX/2 • **MP 
HO  =  HMIN 


-HZD2  =  HO/2.0 

_ H  =  HZDj|_ 

RKT  =  T 
El  =  r 


r»  rvr» 


SUBROUTINE  ASET 

THIS  SUBROUTINE  CALCULATES  THE  COEFFICIENTS  FOR  THE  AXIAL _ 

DER  fvATm“”PACKAGfc  AND 'ALSO  INITIATES  THE  W  ARRAY  AND  THE  W  L 

ARRAY*  _ 

6lMtNSiO^  FC(80>.  NM( 20 ) *  BD ( 10*8 1 , 15 )  *  Adi)  * 

2  »AZ(81) *DD(81) »  PVD ( 8 1 ) .BVD ( 81 ) »  BZD(bl) 

DIMENSION  DVTH( 10*81)  *  OTTHI 10*81 )  »D?TTH( 10.81 )  _ 

common”/ e 7  fc”  nm””/b>  bd 

EQUIVALENCE  (bD(l),  P)*  ( BD< 1621 ) » T )  >  (BD(4861)>  RHO). _ 

1  (  BD ( 729 1 ) »  V).  (BD( 10531) *W)  »  ( BD ( 11341 ) » WZ ) 

2  ♦  (  BD  ( 39 1 1 )  >  DV  TH )  .  (  BD  (  3241J_*  DTfHi*  (BD(4Q51)  .02TTH) _ 

~  EQU I  VALENCfc“TFC  (  33‘)”"TeD  )  ♦  (  Ft  (  39  )  *  St )  *  <  FC  (  33  )  .  DtLV  )  * 

1  C  FC ( 32  3*  XL).  (FC( 34) .  XJ).  C  FC  C 

T~7T27ai1Tf5amT»  (TcTi¥TV  oamT 

3  * ( FC ( 59  )  »  SRD ) >  < FC < 60 ) . SC6 ) » ( FC ( 6 1 ) » SCR ) » ( FCt 62  )  .SIP  )  . 


4  ( FC ( 38 ) *  ZIP)  » ( FC ( 64 ) » BB) * ( NM( 14 ) .  ND ) 
t»<FC<6_3)  *SIP2_>.  (  FC  (  65  )  .  BC )  »  (FC(o6).BZ) 
DO  ”40  T  =  l'.'ND . 


BVD  C  I  )  =  RHO (1*1)*  V< 1.  I  ) 

Tvff(  i r  *  im  »T”*  dvt“h  rn  i  r 

DD( I )  *  DVTH( 1*1 )**2 
BZD (  I  )  =  CJVDf  I  )*6tTh(  1. 1  ) 

IF(RHO( 1.  1  )  )  60.20.20  _ _ 

20""w  (  1  .  I  f~*  (2 •  0  +  S  C6  *  S  Q  R  TF  <  RhOTl”  I  )  |*(V(1.  I  ) **2+DEL2V )  **.  2_50*SRi) ) 
1  /SCR 


Ml 


- x2Tr“-m7rr*vTTrrr«T - 

40  WZ(l.I)  *  W<1.I)*(ZIP  -  T ( 1 . 1 ) ) 

C  INITIALISATION  OF  A-  ARRAY  tMPLOYlNG  W-  ARRAYS 
_  _  A ( 1 )  =  WEDD(RHO)  _ _ ; _ _ 

^wTnt“=''we5'd  TwT 

_ ^ _ cm  =  Bb_  *wint  __  _  _ 

---- (-gam"“”I.6) *“we"ddTp) 

_ A  ( 4 )  »  VZ»  A  (  1  ) _ 

BAB  =  WbDS(AZ) 

BAD  *  WEDS(DD)  _ _ _ _ 


BABE  *  -WEDS  ( BZD )  __  _ _  _ 

'"fTnT~=''bz"*sTp"*b'ad'  .  "  ~ . 

_ FINK2  =  SIP2»BB»BAB _ 

FINK4  *  -(GAM  -  1.0}*BAE 
FINK5  =  BC*WEDD(D2.TTHJ 
A  (  6  j  =  A  (4)  +  WED‘D(DVTH  )’*  SIP’*  BZ 
C ( 2 )  e  BB*WEDD(WZ)  +  FINK  +FINK2  +  FINK4  +  dAd&+FINK:> 

_ A(  8  )  =  A  (  1 ) _ 

C ( 3 )  *  -C ( 1 ) *DELV 

RETURN  _  _  _ _ 

60  CALL  DRAW 
STOP 
END 


C_ _ THIS  FUNCTION  EMPLOYS  WEODLtS  RULE  TO  EVALUATE  THt  INTEGRAL ( 0»2PI 

51  m  EnsTon"  BTeT)  V‘frcT8  oT  »  *nm  ( 2  57 . 

COMMON  /E/  FC  »NM 

£Qu^VAL^Nct  ( fslM  ( 1  S> )  ♦  N1)*(^C!(54)»DTH)  ' 

_ SUM  J  0.0  I 

'  D0"30  ------------  . 

30  SUM  =  SUM  +  38.*B(I)  +75**(b(I+l)  +  b(I+4))  +  30  •  *  (  o  (  ! +<: )  + 

- T - STIV57) - - - 

WEDS  =  S.U*DTn/28fa.*SUM 

RETURN  I 

END 


FUNCTION  WEDD(A) 

this  function  employs  weddles  rule  to  evaluate  the  integral^. 2PI j 

■"  D I  MENS  I  ON" '  A  (Y 6' V§  l")  *""bTa T Y 7“  F'c"("§  oT.'“NM'('2  67 . . . 

COMMON  /E/  FC»NM 


EGU I VALENCt  (NM(15)t  Nl)#(l 
DO  10  I  =  1 »Ni 

-C(54)»DTH) 

10  5(1): 

=  A  (  1  *  I  ) 

SUM  =  0.0 

DO  30 

I  =  1 »N 1 » 5 

• 

30  SUM  = 

SUM  +  38. *B ( I ) 

+  75 •* ( B ( I  + 1 )  +  b ( I +4  )  )  +  50.*(d(I+2)  + 

1  5(1+3)) 

WEDD  =  5.0#DTh/268.#SUM 

« 

RETURN 

• 

END 

IS* 


owonww  mu 

THIS-  ROUTINE  COMPUTES  THE  PARTIAL  DERI VAT IVES  WITH  RESPECT  TO  T 


»  if**1 


(63)  * 


X  ( 4 )  «  BO ( 10*81 • IB ) •  DRHOT ( 10*61 ) *  DTT(10»61)t 
1  DVT (10*81)*  RHO ( 10*81 )  •  DVTH(  10*61 ) *DRH0TH(  10*B1  )*  W410*BD* 


2  DPTH(10»ei)*  DZVTH( 10*81 )♦  D  T  TH  t 10  #e 1  )  •  PI  10*61)*  U2TTM(lC*r 

3  WZ ( 10*81)  •  V ( 10*81 ) •  T ( 10*81 ) »  FC(80>*  NM(20) 


COMMON  /E/  FC*NM  /B/  B0 

EQUIVALENCE  (BD(l)*  P)«  (BD(bll).  OPTH)* 


1  ( BO ( 1621 ) «T ) •  ( BO ( 2431 ) *0TT ) •  (bD(3241).  DTTH ) •  (BD<4031)*  DZTTH) 
2 » ( BO ( 4861 ) *  RHO)*  (B0(5671)*  DRHOT)*  (B0(6481)*  ORHOTH ) • 


3  ( BO ( 7291 ) *  V)*  ( BO ( 8101 ) •  DVT)*  (B0(8911)*  DVTH ) • 

4  ( BD ( 9721 ) ♦  D2V1H)*  (BD< 10531 ) *M) •  ( BO ( 11341 ) «WZ ) 


EQUIVALENCE  ( FC ( 32 ) •  XL).(FC(34).  X J)*  (FC(33)*  RED)* 
1  ( FC ( 41 ) *  VZ).  ( FC(42 ) *  FOAM) •  (FC(36)*  GAM)*  (FC(53 


(64) *68) • (FC(fe3 ) *8C) •  (FC(66)*BZ)*  (FC(62)» 
3  • (NM( 14) *  NO) 


*  VZ*X ( 2 ) 
BE  *  VZ*X(3) 


THE  DERIVATIVE  OR  RHO  WITH  RESPECT  TO  T  -  THE  CONTINUITY  EQUATION 
THE  MOMENTUM  EQ.  -  THE  ENERGY  EQUATION 


DO  40  Is  1*  NO 

DRHOT (1*1)  *  -RHO ( 1 . I ) * ( DVTH (1*1)  ♦  X(l)>  -  V ( 1 ♦ I ) *DRHOTH  ( 1*1 


)  -  BA  +  W( 1.1  )* 

DVT  (1*1)  =  (-(RHOU.I  )*DVTH(l*i  )  +  BB*W  (1*1  )  )  *V(  1 » I  )  - 


2  DPTH(1.I ) /GAM  +BZ#D2VTH ( 1 » I ) ) /RHO (1*1) 

3  DTT ( 1 » I )  =  -V(1*I )*OTTH(l*I)  -BE  +  ( ( i.  -U AM )#P (!•!)*( OVTH ( 1 » 

TriT+'WwTTTTTT—'Br^Trrj? - *"rDVTHTT*T" 

2)**2  +  X(1)*(X(1)  -  DVTHCl*!)))  +  Bb*W(l*I )*SIP  * 


3  V ( 1 » I ) **2 )  /RHO ( 1 » I ) 

RETURN 


END 


_ JJBT"  JiAb'M  (yp  .  Y »  kk  )  *4*ur  '>■•■ 

..JjMENJiON.YPl6i0)A.JtiAJjiiA.FC16fi.[i„tlMJjiClli„WTil0i. _ 

COMMON  /E/  FC.NM 

EQU I  VALENCE  ( FC ( 3 ),H),  (FC(ll).WT),  (FC(lO)t  El).  (FC(9).  tMAX ) 


1  { FC ( 8 ) *  EMIN).  ( NM ( 2 ) »  INDR). 

2  ( NM ( 3 ) .  MALP)  »(NM(20).  NOD) 


( NM { 7 ) »  NSW).  (NM(4)»  MGAM) » 


rA  A.  1' 


*####**#*********#**##**  ******** 

MGAM=-1.0.1  INDICATES  PRcD I L TUR-KUNGE-lOJTTA  OK  CORRtCTuRPnASu 


■mALP-TndTCAY t  V  P E R T “  0 F  RUNG t  KU T T A  PHA S t 
NSW-  PRINT  SWITCH  FOR  ERROR  INDICATION 


,f-  u 


C  El  -  CONTAINS  MAXIMUM  ERROR  FOR  EACH  CYCLE 

C  YP  -  ADDRESS  OF  DERIVATIVE  ARRAY  _ 

y"~-"~ADD”rES"s  c”f””t7e~  ORDI  NATE  ARRAY  . 

Q  ******************************* 


iy  ( MU AKT  40.  10*  60 


C 

C 


THIS  ROUTINE  EMPLOYS  A  FOURTH  ORDtR 
**  **************  *  *****  **  *********  * 

TfTmal"p"'-'"2T‘2‘67Ts”",T5" 

_  <■ =  4  -  MALP 

"do  18  T  ='T,“n'od",io" 

J  =  I  +  K 


RUNGL-K.UTT A  STARTER 


10 

15 


18 


20 

22 


Y ( 1+9)  =  Y( J)  +  H*YP ( I ) 

GO  TO  99 

DO"22”  l"="-r,"NOD7iO" 

Y ( I +9  )  =  Y ( I  +  3 )  +  H/6 • * ( VP ( I  ) +2 • * ( YP ( I ”1 / t  , 

GO”  TO“  59“  . 

*  *******  *  *  *  *  *************** 


.  +2 )  )  +  YP ( I+J ) ) 


"  tf  -i  i. 


C 

C 


ADAMS  BASHFORTH  THIRD  ORDER  PREDICTOR  FORMULA 
*  *  *  *  ************ ***************** 

“DO  4  2~ ” T ~  "  i“»" NO 5T  I”c 

Y(I+9J_  =  Y  (  I  J  +  H/12,0*(23«0*YPl  I  )  -16»0*YP  (  I  +  1) 

"go”  To  99  . . 

*******************  X-**  >***** 


40 

42 


+b • 0*YP ( l+t) ) 


C 


ADAMS  MOULTON  THIRD  ORDER  CORRECTOR  FORMULA 


C 

Jl 

C 

C 


*60*  DO  9~8"~T”  =” i’.“n'65’.To”* ’  ’  * . ”  ' . 

_ Y.i.U_.r__Y_(.I  +  lJ._+_Hy_U!H*_L^_._Q*YF_Li.L..t_liJl?Y_P_U+_iJ__- .YPl.l+ZlJ..... 

E  =  AbSF(Y(I)  -  Y< 1+9 ) )*WT ( RK) 

_ IF ( Y  (  I  )  )  7C  *  80*  7  0  _ 

70  E  -  E/ABSF(Y(I ) ) 

aC._..LE.L.A.-.E_MAXJ._8^j„9Ai..S^ _ 

65  I F ( E  -  EMIN)  98,  87,  87 

*  *  **  *******  *  *  *******  *  *******  *  ****  * 

“rela t"i  ve " er r"or”"c~he~ck.”br~an”ch""t"o_-99  I"ndTc"ate"s~~error~  sm a l l"er*"th"an 

_ ALLOWABLE  ERROR-ADDING  ONE  TO  INDR  INDICATES  VARIABLL  within 

ERROR  ALLOWED 

■a*************************  ******  * 

87'T\DR'VTNDR  +’“l 
I  F ( NSW )  88.  98,  88 

88  "wrTtY  T  6*2  )” 

2  FORMAT  (  1HQ  ,  8HVAR  I  ABLE  .  13,  13HW  I  THIN  'LIMI  TS  ) _ 

GO  TO  98  r 

****************************** 

ONE""MUNDlREl5"7s"~sWfRXcYt5""FOR"”trACTl"V^rA8UE~TAl^Gr^"fiiAN’' Trie 

_ £_R  RJD  R _ L_I_M  J_I  S_ _ 

*V*V*V*V*V*  *****“***  *"*“**“*  *“■****"**"*  ~ 

95 _ INDR  =  INDR  -  100 _ _ 

I F ( NSW  )  97  ,  96,  97  ~ ’ 

97  WRITE!  6,3)  KK  J.3  Si  .  _  „ 

3  FORMAT"?  1H0 ,  8HVArTa”Bl"E  ,  14,  20HQUTSIDE  ERROR  LIMITS? 


i 


o'ri  r» 


*#■»■*#*  *#**#**##*##*  a*#*#########**#* 

£1  CONTAINS  MAXIMUM  tRROR  0CCUR1NG  DURING  Th£  CYCLE 
*"##****V  lViV*V#V#¥*V**V*V#ViV#lT 

98  £1  =  MAX IF ( £  »  £1) _ 

"59  R£TURN 
»j_  END 


DIMtNSION  P(  10*81)  .FC(80)  .  NM  (  20  )  .  BD  (  1 0  ♦  8 1  *  18  )  »  ARK  o00*^) 

. common-art' .  .  . 

COMMON  /E/FC»NM/B/BD 

Equivalence  (fc(d,t)»  (bd(1>»  pi.  inm(13)*  iuj  " 

1__»_  (NM(14),  ND)*  ( NM  (  1 5  )  *N  1 ) 

_X M  I  =  XMA 

sum '=  xm'i . 

_ DO  10  I  =  2«N1 _ 

SUM  =  SUM  +  P ( 1 »  1  ) 

I  F  (  P  ( 1  *  I  )  -  XMA)_  5  1 1 0  »  3  _ 

=“p-riViT 

GO  T_0_  10 

5  '  '  'TFipTYrrn'-  xmIt  7V~y67T6 

__7 _ XMI  =  PjJjJJ _ 

10  CONTINUE 

XN  =  N1  _  __ 

SUM  =  SUM’/XN 

ART  (  I U » 2  )  =  (XMA  -  XMI) /SUM  _  _____  __ 

art  ( iu~,  rr  *=  y~ .  . . . 

NN  =  Nl/4  +  1 _ 

ART( I U  »  3 )  =  P ( 1 » NN ) 

RETURN  _  ■ _  _ 

END 


*?rv  >»* 


C 

c 


SUBROUTINE . TKPREIT  :  ~  ~  " 

EMPLOYING  THE  RESULTS  OF  THE  INTEGRATION  THIS  ROUTINE  ATTEMPTS 

comPuTe”! H^E'  f he t”a"""d"erTv^tTvTs~”6P*  f Pit  n+I  "“annulus" 

DIMENSION  FC  ( 80 ) »  NM(20J.  faD (10 *61 » lb ) »  DRHOTH ( 10*81 ) » 

D2vTH( 10*61) .  DTTH(lo»81). 


TO 


T  dvthI lo.bl ) * 

2  DPTH ( 10 *6 1 ) » 

T’t  no.  §TT aa"9  oTV  ab  ' ("9  57 » 

COMMON  /E/  FC »NM  /B/BD 

E'5uTv7fLlTic"r"(BDn“*""p37"TiD(8TlT*“b’PTTTj"rTB'D(T67r"r»"' 

1  ( BO ( 3241 ) *  DTTH ) .  (B0(4051).  D2TTH),  IBD(4861).  RhO). 


02TTH ( 10»6l ) » 

RHO (10*61)*  V ( 10 *61 )  » 


T) 


P( 10*61) * 
"AC"(90") 


UVTH) 


) 


3  (  E.D  (  6481 )  »  DRHOTH).  (BD(72Y1),  V).  (60(8911). 

3_(BD(9721).  D2VTH).  (FC(84)»  DTH )  »  (NM(4)_»^  MGA 

4  "*TNMrri4l7ND )  7"Tnm  (Ts ' n  IT^TfcT  ' ToTi  bTTT"  (  fc  77  i  f*  dsu7""'~ 

5  . (NM( 16) *NZ) ♦ (NM( 17) »N2 ) »  (NM(18)*N3)»  (NM( 19) »N4) * (NM( 12) 


N5  ) 


10 


IF 

J  =  10 


(MGAM)  10.  10.  12 


12 


13 


K  =  1 
GO  TO 
J=  T" 
K  =  2 


13 


AA(  1  ) 
AA  ( 2  ) 


RHO (  J  »  N5 ) 
RHO ( J  »N4 ) 


A  A  (  3  )  = 
AB  ( 1 )  = 
'ab  ~(Tf“=" 
AB  (  3  )  = 

ttcttts" 

AC ( 2 )  = 


RHO ( J  » N1 ) 
V( J.NS  ) 

,7("J»T\74") 

V(J,_Nl_) _ 

7T”J»7\”5~) 

T  (  J.N4) 


AC ( 3 )  =  T(J,N1) 

DO  20  1=  4.NZ 


AB ( N3 )  =  V ( J  *  3 ) 

.ACJ_N2_)___= _ IJ.J_._2J_ 

AC ( N3 )  =  T ( J  *  3  ) 


j 


*r- 

>  •* 
**  ■ 


C  FOURTH  ORDER  STIRLING  CENTRAL  DIFFERENCE  FORMULA 

£_  *******************  *#j*  ******  #***•»■#**#'•#■#*_#* 

™--~  ----------  - 

_ _DRH_0_T_H_(_J_._I_J__f__i_AAj_I_+lJ_-a_.  p_ti_AA.(_I_+2.)._— _AAJ.l+4JJ__r__A_AJJ.+_3_J-)_/p_^ _ 

~  dvThTjViT"  ="Ta  bTI"+  l)  ’  - "  q~~  o’*1"a  5Ti"+ 2T-a’b"u  ' +  41"  )"-"a  d~u  +" 3 17/15  2’ 

DTTH(J.I)  =  (  AC  (  1  +  11-8. Q*(AC(  I +2 ) -AC (  I  +4 )  )  -  AC(I+5))/D2 _ 

D2VTH ( J  ,  I )  =  ( -AB ( I  +  1 ) +  16.0*( AB ( 1+2)  +  Ab(I+4))  -  30.*AB(i+3) 

_ -L__ .rA3XI±5J_iyH^l _ 

25  D2TTH(J,I)  =  (-AC(I+1)+  16 . 0* ( AC ( I +2 )  +  AC(I+4))  -  30.*AC(i+3) 

_ J^LASJLl±5J.l/CSfl _ 

35  DRhCTH(J.ND)  =  DRHOTH(J.l) 

_ DVTH  (J.ND)  =  DVTH  (J,l) _ 

D2VTH  (J.ND)  =  D2 VTH ( J . 1 ) 

DTJH  (J.ND)  =  DTTh(J.l)  • 

D2TTH( J.ND)  =  D2 T  T  H  <  J » 1 ) 

_  _ _  DO  50  JL  =  1»ND  _ _ _  ■ _ 

P(  J.l  )  =  RHO (3,1  )  *  T ( J  ♦  I ) 

50 _ DPTH(J.I)  =  RHO(  J.l  )»DTTH(  J.l  )  +  T(J*I)»  DKhOTn(j.l) 

RETURN 

_ IHOL _ JL41 _ _ 


- SJBffOUTOE^HTrr - 

DIMENSION  FC ( 80  )  , ' NM ( 20 ) .  BDU21P0)  _  . _ ' 

(T6^OT"7r7  rrr'NM  7b?  sir  "*  "  : 

EQUIVALENCE  (NM ( 5 ) *MPTN ) * ( NM( 4) »MGAM) •  <NM<3) ,MALP) . 

1  {nM(I).IM).  <P£(3)»H),  <FC(7) »H262) »  ( NM { 10 ) »NV ) ,  (FC(4).H0)» 

2  IFCI1).T)»  ( NM ( 2  )  *  INDR).  (FC(10).  El).  ( NM ( 8 ) .NCOU ) » 

.  3  ( FC ( 6)  .  HMA'X)  »*7nm7"6~)  »MP“fS)”.  ( FC~f57*HMTN )  .  (PC ( 9T.’  cMAT)  " 

4  t  FC ( 2  > . RKT ) .  ( NM ( 20 )  .  NOD ) 

“  iViViT#*#ViV###*"*V**Vi#*'*iV*Vi#####V»V*#*'*****i**###***##*****#*# 

C  SHIFT  UPDATES  THE  VARIApLES  SO  THAT  PROPtR  POSITIONING _ ;_ 

C  OCCURS  DURING  THE  PROPER  CYCLE 

C  I.***##*###**#*#****#**#***#**#*##***#*#  a-*##*######*##*#*#***###*# 

IF  (  MPTN)  15".  I0“.“l5 

_ 10 _ MPTN  =  MPTS  _  '  _  __  _ 

15  MGAM~"=  -MGAM  ~~  . 

I F (  MGAM  )  80.  20.  60 
20  I F (  MALP  -  2)  25.  48.  4b 

_ _ _ 25  IM  =  1M  +  1 _ _ _  ___  ____  _  _ 

MPTN  =  MPTN  -  1 

_  _ NCOU  =  NCOU  +.1 _  _ 

'  .  I  FI  3  -  IM)  35.  30.  30~  "  ’  . . 

_ 30  MALP  =  4 _ 

H  =  HZD2 

_ GO  TO  40 _ _ _ _ _  _ _ _ 

--- 

40  DO  42  J  =  1.15 

#  g-j-Q  +  j  .  - 

_ NZZ  =  NZ  +  NOD  -  1 _ 

DO  42  I  =  NZ .NZZ .10 

_  _  . _  _  _  BDU  +  1)  =  BD(.I+3J _ _ _  _ 

BD(l  +  2)  =  BD ( 1+4 ) 

_  _  _  BD(  1  +  3  L  *  tJD(_I+5 )  _ _ _ 

BDTi+4)  =  60(1+6)  .  7 

_ BDU+5)  =  BD  (  I  +7  ) _ 

42  BD  C I  )  =  BD ( I +9 ) 

_ REHASH _ 

Q  *******#***#*#********#*#•**#*###*###*#*■«.**#*********##***#*#*#*** 

_ £. _ SJlI  F.IJLN.G._0_CpjRS-AF-XER-PR£DXC10  P_X_YCLE^1  ACQMP  LE  IE _ 


C  *#***##  ****************  ****#+**#*#*#*•**#*#***#**##*****#* ******** 

_ 46  MALP  =  MALP  -  1 _ 

IF;  MALP  -  2)  55.  50,  55 

_ iQ..tL-r-jdfl. _ _ _ _ _ 


GO  TO  65 

_ 

GO  TO  65 

60  T  =  T  +  H _ 

65  DO  75  J  =  1,15  '  ' 

_ CJZ__=_-LJr_lJ.-*-.e-LO._±„l _ ; 

NZZ  =  NZ  +  NOD  -  1 

_ D^_._2A_JL_=-£LL»A2jLi-L0. _ _ _ 

BD( 1+8)  =  BD ( 1+7 ) 


C  THE  FOLLOWING  STATEMENTS  EXIST  FOR  CHANGING  THE  STEP  SIZE 


Q  *##*####*#**#**###*******################*##############*#*#*###* 

6 O'TfTi'nDRT' l'7b  ,"l  35 V 125 . . . 

125  MPT N  =  MPTN  -  1 


130 

IM  =  IM  +  1 

INDR  =  0 

El  -  0.0 

NCOU  =  NCOU  +  1 

-  - - 

. 

CL 

RETURN 

***************************************************************** 

DOUBLE  THE  INTERVAL  IF  POSSIBLE 
************************************** 

135 

140 

IF( IM  -  6)  125 ,  140.  140 

KZ  =  MPTN/ 2 

C 

I r ( 2  +  K.Z  -  MPTN)  150  .  125  .  150 

*************************************** 

c 

c 

IF  MPTN  IS  ODD  THERE  IS  AN  EVEN  NUMBER  OF 
************************************ 

POINTS  LEFT  TO  COMPETE 

150 

155 

IF ( (H-HMAX ) /HMAX )  15b,  125.  125 

I F ( MPTS  -  1)  125,  125,  160 

160 

MPTS  =  MPTS/ 2 

MPTN  =  MPTN/2 

IM  =  5 

DO  165 J  =  1,15 

NZ  =  (J-l)  *  810  +  1 
NZZ  =  NZ  +  NOO  -  1 


DO  165  I  =  l\!Z  *  NZZ  *10 
oD  (  I  +  1 )  =  pP  (  1+2) 
dD(I+2)  -  PD  (  I +4  i 


_ _B0_(_I+3J. 

165  BD  (  I  +4  ) 
H  = 

GO  TO  130 


p_DJ_I_+_6  J 
BD<  1+8  ) 


C  HAlF  THE  sTc.P  SIZE  IF  PuSS i cLp 

Q  *  ■*  #  &  &  &  ie  St  x  -fc  ir  •&  *Jc  Hr  Si*  ir  *  x  *  ic  ->  ye  Hr  •*  v  # 

111  1 >"  H~ h -~H  M I  N  )  / h M I  nT ’ "i 25 Y'l'ioT 16  o’"" 

_io0  N  f  LCGF_(  E  1 /EMAX  +  1  .  )  /  ._69  3  147 lo 05 

’""55”  185  I  ”=  1 » N 

_ K  =  I _ 

IF  (  H  /  2  •  **K.  -  HMIN)  195  .  190  .  185 
1S3  CONTINUE 

"~"Tso”’’"n"  =  "T" 

GO  TO  200 
195  N'VT-I 

2C0  MALP  =  4 _ 

MG AM  =  0 

I  F  ( _ I M  -  4)  210.  215,  210 

‘2To“""t'V"t""h  . 

_ R<T  =  T _ 

NTT"::'  T" 

GO  TO  216 
TT5  ;  =  RK.T 


MPTN  =  MPTN  +  4 


NN  =  5 

2a8_  MPTS=MPTS#2#*N 
MiP  T  N =MP  T~N#~2#"*N~ 

DO  222  J  =  1,15 _ 

NZ  =  (J-l)  *  810  +  1 
NZZ  =  NZ  +  NOD  -  1 
DO  ~22  2”  T"=”  nZ  » NZZ  *1  6” 


4 


Kr  *  I  +  NN 

222  ^D(I  )  =  BOCK) 
nO/  =  "'h7  2".  **N_ 
H/D2  =  HO/2. 

H  =  hiJ2 

Am  =  o 
rTNDR“""="b" 


J_  RETURN 
‘END 


SUBROUTINE  DRAW 

DIMENSION  ART  (  800  ♦  3  )  ♦  Ad  <  82  .40_)_»  _f  CJ_bOJ. *NM (  ^0  )  »_bC  dpj_3.).» . . 

l  "acTe  if" 

CUMMoN  AKT  *  ao 


COMMON  /t/  rC,  NM 

tOO  1  VALtNCt  ( NM l 1 i ) »  I  M )  .  (  NM  (  10_)_«NN  )___»_  J  NMJJ  4_)_  •  ND_)_  >  I 

"f YP C* ’ I "n ft" Gt K " ’  oC V  00  . . . . . 

J  L  -  Ii*i 


.  =  N 1 

X  =  6.2832/XN1 


AC (  1  )  =  1.000 
DO  4 1  JM  =  2. NO 


42  AC(JM)  =  AC(OM-i)  +  IX 

NZZ  =  ARTUZ,  1)  +  1.0 


MAXIMUM  AdSC 1 SSA  /ALUfc 
XMA  =  NZZ  £  VM  =  C.O 


DO  20  1  =  1,  JZ 

AR  T (  1*1)  =  7.0/XMA*ART ( I .1 )  +  1.0 


IK (ART! 1 .2)  -  VM)  20.  20*  10 

10  YM  =  ART ( I ,2 ) 


20  CONTINUt 

YM  =  10.0*YM 


NZZ  =  YM 
YM  =  NZZ  +  1 


YM  =  YM/10.U 
YM I  =  0,0 


flC (  1 )  =on  T  IN  RA£dC(2)=  fanDlANS  £  oC(3)=  bH 
sO  (  1  )  =orl  (  mMaX-PMS  dl>  (  <£  )  =oh  1 N  ) /PAVu  »  do(j)  =  otic 


^aLl  ukMKn ( mkT ( 1 » 1 ) , aR7 ( 1 » 2 ) » XMa * 0 • u  » YM » YM 1  * JZ » oC » au «  -l) 

YM  =  l.o  x  YM/2.0 


YM I  =  2.0  -  YM  x. 0000001b 
oO  =  C • COu 


DO  232  I  =  l.JZ 

IF  ( ART ( 1 ; 3 )  -  UP)  232,  232.  231 

231 

232 

UO  =  ART< I .3) 
CONTINUt 

I  F ( UO  -  YM)  234, 

234,  ^33 

23  3 

NZZZ  =  10.*UO  + 

1.0 

UO  =  NZZZ 

YM  =  00/10.0 

234 

CONTINUE 

8D( 1 ) = >OMAX IMUM  $ 

bD(_2j  =dHPKESSUKE  £  bD  (  3  J  =  bn  NOPE 

CALL  GRAPHIARTI 1, 
OZ  =  ND 

1 ) ,ART< 1,3) ,XMA,0.0»YM,YMI , JZ , oC . dD , - 1 ) 

oC  (  1  )  =  fcnTHtTA  IN  s>  dC  (  2  )  =  8H  RADIANS  £  oC  (  3  ) 
bD(l)  =  bnPRtSSURt  £  dD(2)  =bH  WAVE  PE  £  BD(3) 


DO  62  L  =  1 » NN *  6 
NZN  = 


CALL  GRAPH  ( AC ♦ AB ( 1 , NZN ) . 7 . 0000 » 0 • 0  * YM , YM 1 »  JZ»bC,6D»+l) 
DO  60  I  =  1.3 


JL  =  I  +  1 
NJN  =  NZN  +  I 


60  CalL  uKAPn  ( mC . Ab ( 1 » N ON ) . 7 . 000 0 » 0 • 0 » YM . YM I . — JZ * dC . dD * J L ) 
NZN  =  NZN  +  4 

o 2  ‘call " gr a'pFTTacVab ' 7T 7n, zn  i TT.'o o"667o’.‘67y ymT» jz"»" d c .' ao 7- 5”j ' 

_  RETURN 

'END'" . .  ~  . .  '  ' 


SUBROUTINE  GRAPH! AB»CR,XM»XMI #YM#YN' tJLt  »«.** 

DIMENSION  ABI800) *  OR(800)i  BCX (3 ) *  BCY ( 3 ) •  FC ( 80) »NM ( 20) •  6(400) 

TBBW0JT‘7E7Fr^M - 

EQUIVALENCE  ( FC ( 32 ) . XL )  •  (FC(33)»  RED)*  (FC(34)*XJ)» 


f  worn  lOMUMl—flai  Hl»i:n 


2*  ( NM ( 10 ) »N J ) 

DatTYmM) . 

UNIT  17  DESIGNATED  FOR  PLOTTED  OUTPUT  ( EQUI P , 17«** ,SV ; 

- irmrrTrTr-9- - - - 

8  CALL  PL0TS(  G.  400,  17) 


M  Z  = 

MZNM  =  NJ  -  1 

SCALE  6  R  fi  IN A T E" VALUES  0F"InC0MI  NG  DATA 
9  NPT  =  XABSF.JZ) 


(MS) 

DIF  =  9.0/ ( YM  -  YMI ) 


DO  14  I  =  1,  NPT 

14  OR  (  I  )  =  (OR!  I  )  _r  YMI  )*DIF  +  .5000 

DATA  SC ALt5  FOR"  ~~N I_N E’"l NC h“" O RD'lNA T E  "FRAME 
IF(JZ)  70,  15,  15 


15  NZ  = 
ROW  =  2 


SIGMA  =  ( YM  -  YMI ) 

SIG  =  1.0  4  SIG1  =  .1 


I F ( SIGMA  -  SIG)  21,  21,  20 
YB  =  2.0*SIG1 
"GC“fO~  25" 

21  IF(SIGMA  -  • 50*S1 G )  23,  22,  22 


22  Yo  =  SIG1 
GO  TO  25 


=  NZ  +  1 

G  =  J  •  0  /I Q.  Q.**N  Z.  _  _ _ 

gi  ="sTg/i"o".o  . .  . " . .  . 


25  YZ  =  DIF*YB 

L _ S_CAi.ljiG__PA_RAM£IERS__HAV£:,d£EN..SET^FQR^ORJlIIYAT£ _ 

C  SET  PARAMETERS  FOR  ABSCISSA  SCALING 

_ XLA _ f _ AMI. _ 

XAB  =  1.0 

-  XM 


STG  =10.0  4  SIGMA  =  XM  -  XMI  4  ZM  =  1.0 

_ 21. _ _ iQj._3LQj-_2J. _ 

28  ZM  =  ZM  v  1.0 

_ .QSL.J.Q.21 _ _ _ 

30  XZ  =  D1F*ZM 


C  BEGIN  ABSCISSA  PLOT 

_ CALL-MUMHER—L-Ai— Q-^3.5-i--iD2j  XLAt  Q  tfl.i_AtlEfu.lJL _ 

CALL  PLOT  (XAB,  0.5,  3) 

_ ,vALL_PJLQI__USABjL_iL*£i_^ _ _ _ 

XAB  =  XAB  +  XZ 
XLA  =  XLA  +  ZM  4  MM  = 


I F ( XAB  -  8.0001)  34,  34,  35 

_  34  CALL  PLOT(XAS,  9.5,  3)  _ _ _ _ _  _ 

CALL  P~L0T  (  XAB ,  C .17" 2*7  . 

_ _  X  =  XAB-.07  _  _ _ 

CALL  NUMBER!  X,  ^ 3"5",‘^0‘7‘,"x"la7oVo ,  4HF4.1) 

XAB  =  XAB  +  XZ  S  XLA  =  XLA  +ZM  4  MM  =  1 
IF  (XAB  -6.0001)  31,  31,  35 

_ 35  X  =  8.0  +  XZ  -  XAB _ _ _ _ _ 

ZETA  =  XZ/2.0  $  ZEB  =  "x AB  "7 ETA  SXLA  =  x‘la"-Z”M/_2  ."6+7oo5o~oT 


IU1 


Hi 


~bJ3K Ci'JTWLXAPh  [a b ."OK , x M , xTQ  , Y M *Tm f»  JZ ,  B Z X  »“3T? »  nS ) 

DIMENSION  Ab  (  300  )  *  OR  (  800  )  *  BCX  (  3  )  *  BCY(3)»  FC  (  80  )  »NM  (  20  ) »  G(400) 

— loro— 
0020 

CuMMvN  /E/r C  »NM 

jC 

ECU i VALENCE  ( FC ( 32 ) *XL )  ♦  ( F  C ( 3  3  >  »  RED).  (FC(34).XJ). 

0040 

rnrrcrb 6 ) »  o  am  ) .  o z rsv  yvuti' v ) .  tftt iwn'3")"."'i^i 

oo^o 

(  .. 

2.  (\M(10)»NJ) 

DATA  (M Z=C) 

6  0 

UNIT  17  DESIGNATED  FOR  PLOTTED  GuTPUT  ( EOUI P , 17=** ,SV ) 

0070 

1  r  (  MZ  )  9,  8  »  9 

cu 

• 

0 

CALL  P LOTS (  G»  400.  17) 

VO 

MZ  =  1 

i  00 

MZNM  =  NJ  -  1 

SCALE  ORDINATE  VALUES  OF  InCUm  I  No  DATA 

01 1  u 

9 

NPT  =  XAbSF(JZ) 

U  u 

NOT  -  XAbSF(MS) 

1O0 

DIF  =  9.0/ (YM  -  YMI ) 

140 

DC  14  I  =  1,  NPT 

0130 

14 

OR  (  I  )  =  ( CR ( I )  -  YMI)*D1F  +  .5000 

01  bO 

DATA  SCALlD  FOR  NINE  INCH  ORDINATE  FRAME 

0 1  70 

IF(JZ)  70.  15.  15 

17b 

15 

NZ  =  0 

ROw  =  2.5 

lot 

19U 

0200 

19 

I F ( SIGMA  -  SIG)  21.  21.  20 

210 

20 

YB  =  2 • G*S I G 1 

220 

GO  TO  25 

230 

21 

IF (SIGMA  -  • 5C*S I G )  23.  22,  22 

0240 

22 

YB  =  SIG1 

2b0 

* 

GO  TO  25 

2b0 

/ 

NZ  =  NZ  +  i 

2  70 

l 

SIG  =  I.u/10.0**NZ 

200 

SIG1  =  SIG/IC.O" 

290 

GO  TO  19 

300 

•  25 

YZ  =  DIF*Y5 

3  iC 

“ 

SCALING  PARAMETERS  nAVE  BEEN  SET  FwR  ORDINATE 

0320 

SET  PARAMETERS  FOR  ABSCISSA  SCALING 

03b0 

J(LA  =  AMI 

0340 

XAB  =  1.0 

0350 

DIF  *  7.0/iXM  -  XMI ) 

0360 

SIG  =10.0  $  SIGMA  =  XM  -  XMI  i  ZM  =  1.0 

0370 

_ 21. 

I F (  SIGMA  -  10.0*ZM)  30,  30,  28 

j23_o_0  „  . 

28 

ZM  =  ZM  +  1.0 

390 

GO  TO  27 

_ 4_QG _ 

30 

XZ  =  D I F*ZM 

410 

31 

X  =  3  -.07 

_ 420 _ 

3EGIN  ABSCISSA  PLOT 

_ CALL  NUMBER  (  X,  0.35.  .07,  XLA »  0 .  Q ,  4HF4.1) _  _ 

CALL  PLOT  (XAB,  0.5,  3) 

CALL  PLOT  (XAB,  9.5,  2) 

0460 

XAB  =  XAB  +  XZ 

470 

XLA  =  XLA  +  ZM  3  MM  =  2 

0400 

IFvXAo  -  6.0001)  34,  34,  35 

0490 

*  34 

CA^l.  ?LCT(XAd*  9.5,  3) 

C500 

CAlu  PLOKXAB.  0.5,  2) 

0510 

X  =  X Ad— • 0  7 

520 

Ill 

CALl  NuMoER ( X ,  .35, .07,  XLA.O.Q,  4HF4.1) 

0330 

121 

XAB  =  XAo  +  XZ  S  XLA  =  XLA  +ZM  $  MM  =  1 

0540 

Ir  (XAB  -6.0001)  31.  31,  35 

ObbO 

35 

X  =  8.0  +  XZ  -  XAB 

560 

ZETA  =  XZ/2.0  $  ZEB  =  XAB  -  ZETA  SXLA  =  XlA-ZM/ 2 .0+. 000001 

_ 147 _  _  _ 

05  70 

IF (Z£b  -  8*0001 )  37,  37*  38 

XAB  =  zea 


GO  TO  (31.  34)  MM 
ZEd  =  Z£B  -  XZ/ 6.0  SXL  A  =XI 


(  ZEB  -  8.0001)  37,  37,  - 
WRITE  OUT  PARAMETERS  PERT  IN; 


CALL  SYMBOL (  6.3,  9.0,  .07- 
CALL  NUM6ER<6.875,9.0»  .07, 


CALL  SYMBOL (  6.5,  8.73,  .0' 
CALL  NUKdER(  6.873,8.75,  .( 


CALL  SYMBOL (  6.35,  8.5,  .0' 
CALL  NUMbER (  6.875,8.5,  .0' 


CALL  SYMBOL (  6.4,  8.25,  .0 

_ CALk.  NUMBER  (  _6_.  87 5_,___8_.  2 3 ( 

"CATL'SYMBOL  T  6*. 4~3  »  8  V  0_0  * 
CALL  NUMbER (  6.875,  6.0, 


SXL  A  =X  L A  -ZM/6. 


37,  37,  40 

PERTINENT  TO _ Th 

9  *"dr "  i  o"  7  v '  3hU~=" r 

.0,  .07,  XL,  0*0, 

8.73  ,  . 07 ,  3  h  J  - , 
,8.75,  .07,  XJ,  0 


8.5,  .07,  6HDELV 
,8.5,  • 0  *  DELV ♦ 


8.25,  • 0  7 ,  3HRED 
,  8.23, .07,  REO, 


.07,  4H6C 
•  0  /  ,  SC, 


0  +.00000015 


IS  SOLUTION 


0,  3) 

4HF6.3) _ 

0_*  3*) . 

.  5HF11.9) 

=  ,  0 ,  6  ) 

0.0*  4rtF 3.2) 


=,  0.0,  3) 

0.0, _ _4H_F_6._2.J_ 


>»  0.0,  4) 


CALL 

CALL 

SYMBOL ( 
NUMBER ( 

6> 4,7.73,  *07, 
6.875,7.7b,  .07 

BnGAM  =,  0*0*  3) 

,  GAM,  0*0,  4HF6. 

2J_ _  _  _  _ 

CALL 

SYMbOL ( 

3.0,  0*1*  *14, 

BCX  ,  0.0,  24) 

abscissa 

HAS  NOW 

BEEN  LABELED 

,  BEGIN  ORDINATE 

LABEL 

48  YLA  = 

50  MM  =  3 

:  YMI  S 

YAB  =  0.5  S 

MM  =  3 

0710 

AUk _ 

0730 

0740 


CALL  PLOT  (  8.0,  YAd »  3)  t  CALL  PLOT  (1.0»  YAb,  2) 

Y  =  YAB  -.035  __  _ _ _  _ 

TfTnz  i  T  'rd3",”T64TTo"5 

.07,  YLA ,  0.0,  4HF8.1)  * 


CALL  NuMdER (0.7- 


34 

GO  TO  108 

CALL  NUMBER ( 0.7 ,Y,  .07,  YLA,  0.0, 

4hF5. 2  ) 

84 

083 

GO  TO  106 

86 

35 

IFCNZ  -  2)  106,  106,  107 

087 

36 

CALL  NUMBER (0.7,  Y,  .07,  YLA,  0.0, 

4HF6.3) 

w  - 

> 

GO  TO  108 

o* 

37 

CALL  NUMBER (0.7,  Y,  .07,  YLA,  0.0, 

4HF7.4) 

09"0 

)_8_ 

MM  =  MM  -  1 

_ il_ 

I F ( MM  -  1)  58,  56,  56 

>6  YAB  =  YAB  +  YZ _ 

IF (YAB  -  9.50001 )  109,  109,  60 

1SL _ Y__r__YAB._r.#_Q3J2. _ _ 

GO  TO  54 

>8  CALL  P LOT (1.0,  YAB,  3)  S  CALL  PLOT(6.0»  YA 


+YZ  S  YLA  =  YLA  +  YB 

B__-__9._5P_Q.0JL1 _ _P_Q_»__l»_Qj__6_Q _ 

YAb  -  YZ/2.0 

-  9.50 00 1_)_  62  ,__62_»  65  _ 

» 2  YLA  =  YLA  -  YB/2.C 

GO  TO  50  _ 


LABEL  ORDINATE 

.5  CALL  SYMBOL (0.6, 3.0,  .14,  BCY »  90.0,  24) 


Nw,  .  PLOT  CURVES 

0  CALL  °LOT  (  AB(1)  ,  OR ( 1 ) ,  3) 


DO  72  I  =  1,  NPT 

2  CALL  SYMBOL  (AB(I)t  OR ( I )  ,  .05,  NCT ,  0.0*  -2) 
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„2A4l. 

0920 

_0  S.2&. 

940 

0930 


0960 

_Q9_7p_. 

960  ' 
_09_9p__ 
"looo- 
0 


0 

IF  (  MS  )  8U,  82,  82 

call  plot< io.o,  o.o,  -; i 

1090  Or 

2 

RETURN 

ROW  =  ROw  -  *25 

MZNM  =  MZNM  +  1 

CALL  SYMBOL (7.25,  ROW, .05,  NCT,  0.0,~1) 

A 

CALL  NUMBtR (  7.5,  ROW,  .07,  OR(82),  0.0,4HF7.3) 
RETURN 


END 
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